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ABSTRACT 



A coupled ice-ocean numerical model is developed which 
improves the simulation of the annual cycle and interannual 
variations in ice cover in the Arctic. The model is a further 
development of the work by Semtner (1987) . Although the 
accuracy of the simulated ice concentration is increased, the 
annual cycle of ice coverage is still exaggerated. Several 
experiments are conducted to determine the importance of 
incorporating a fully interactive ocean, to select an optimum 
strength parameter for use in the ice rheology, to investigate 
the model's sensitivity to changes in the albedo of the frozen 
surface and to determine the relative importance of the 
various dynamic and thermodynamic forcing mechanisms. The 
regional dependence of these mechanisms and an assessment of 
two statistical analysis techniques used to measure model 
improvement are also examined. 

Inclusion of a fully prognostic ocean component vice a 
ten-year mean ocean cycle in the model improves the 
correlation of simulated ice concentration fields with 
observed data. This is the case for all regions in the 
Arctic; for both the annual cycle and interannual variations 
of the ice cover. A reduced strength parameter value, 
P*=hxl0 4 , is found to improve the simulation of the ice 
thickness distribution with increased overall thickness and 



L— 



* * , v. 



better compression north of the Canadian Archipelago and 
Greenland. In contrast to results using ice models without 
a fully prognostic ocean component, this model is quite 
insensitive to changes in the frozen surface albedo. 
Exceptions are evident where the ocean heat flux into the 
mixed layer is small and the ice is thin. 

At the spatial (110 km) and temporal (monthly) scales used 
here, the heat provided by the ocean appears to be the 
dominant mechanism controlling the position of the ice edge 
and the extent of the ice pack. Within the pack, it is the 
dynamic forcing and, in particular, the wind forcing which 
controls the ice thickness and thickness distribution. The 
ocean circulation below the mixed layer appears to position 
the heat underneath the MIZ. The MIZ is also the region where 
the ice thickness tends to decrease through divergence. The 
linkage between the subsurface heat and the thinned ice cover 
is apparently controlled by conditions at the surface and the 
resulting response of the mixed layer. 
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INTRODUCTION 



Interest in the complex geophysical processes occurring 
in, above and around the Arctic Ocean continues to increase 
as the economic, strategic and political values of this region 
become more important. In the case of economic concerns, vast 
oil reserves are being exploited in this area and the search 
for more reserves continues. Much of Western Europe, northern 
Canada and the USSR are supplied via shipping routes which are 
dependent on weather and ice conditions in the Arctic regions. 
Arctic waters also contain some of the most valuable fishing 
grounds in the world; however our ability to use them is again 
dependent on the state of ice, water and weather. 

The Arctic Ocean spans the shortest distance between North 
America and the USSR and is therefore of obvious strategic 
importance. However the extreme environmental conditions 
encountered there make military operations very difficult and 
expensive. Naval vessels used in northern waters must be 
properly designed if they are to function with any 
effectiveness and their performance will normally be 
significantly degraded. The requirement for an ice breaking 
capability is obvious. In addition, the highly variable 
temperatures, surfaces, and material compositions of water, 
ice and snow have a dramatic effect on both the air and water 
boundary layers. This in turn affects electro-magnetic, 
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electro-optical, and acoustic transmission and reception. 
Sea-ice is a strong source of ambient acoustic noise and radar 
scatter. Detection of surface or low level air targets by 
radar or submarines by passive or active acoustic means 
becomes very difficult. 

Several countries border the Arctic mediterranean and have 
strong political motives for establishing or extending their 
sovereignty and/or economic control as far offshore as 
possible. The presence of several "chokepoints" within the 
region provides the opportunity to influence operations in the 
entire Arctic Basin by controlling only limited areas. 

Obviously, a comprehensive knowledge of the Arctic 
environment, including the ocean, atmosphere and sea-ice, and 
a capability to predict this environment with some accuracy 
would be of great value. The expense of establishing a 
complete observational network is prohibitive. An important 
role must be played by numerical modelling which has the 
potential to provide both a full analysis of the current 
conditions based on limited data, and a prediction capability 
for future conditions. 

A. NUMERICAL MODELLING PERSPECTIVE 

It is recognized that the atmosphere, ice and ocean 
interact in a highly complex nonlinear fashion. Consideration 
of one of these three components cannot be realistically done 
in total isolation of the other two. However a complete 
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understanding of all the physics involved in describing each 
parameter and its interaction with the others is lacking. 
Therefore our ability to mathematically describe the physics 
involved is also lacking. Furthermore, when significant 
theoretical advances have been made in the past, memory and 
speed limitations in computing facilities have been an 
obstacle. Even with today's supercomputers, compromises must 
be made to balance computing requirements (memory, expense and 
time) against the accuracy desired. In order to minimize 
these compromises, the three component problem has generally 
been broken up into its three individual components of 
atmosphere, ocean and ice. The two components which are not 
being specifically examined in a model are then parameterized 
in some fashion or they are represented by prescribed 
empirical data. 

Numerical modelling was first applied to the atmosphere. 
Impetus was provided by public as well as military interests. 
Great strides have been made by meteorologists in predicting 
the weather by numerical models. Numerous atmospheric 
numerical models of varying complexity have been developed and 
are currently in use by the various National Weather Centers. 
However, despite a great deal of effort and expense, 
deterministic prediction of daily weather fluctuations is 
presently limited to five to ten days. Theoretical arguments 
suggest that the ultimate practical range of such predictions 
is about two weeks. 
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Numerical modelling of the global oceans requires a 
thousandfold increase in computer power. Compared to 
numerical modelling of the global atmosphere, the horizontal 
resolution of the global oceans must be ten times as great, 
and the numerical time step therefore ten times as small, to 
accommodate the smaller radius of deformations and faster 
gravity wave speeds. Ocean modelling is further complicated 
by irregular borders, islands, widely varying density 
structures, internal and boundary mixing processes, complex 
and varied boundary layers and a marked paucity of 
observations. Nevertheless with the major contributions of 
Bryan (1969) and numerous others such as Takano (1974), 
Semtner (1974), Cox (1984) and Semtner and Chervin (1988), 
viable three-dimensional, baroclinic numerical models have 
been developed which simulate the state of the world's oceans 
quite well. These models are currently in wide use in 
circulation studies, tracer distribution analysis, etc. 

Numerical modelling of the ice has also progressed 
significantly over the past twenty years. The physics of this 
problem are now considered to be fairly well understood. 
Prediction of the ice both in areal extent and thickness 
distribution, is known to be a function of dynamic forcing 
(both internal and external) and thermodynamic forcing. The 
work of Hibler (1979, 1980) has provided details of the 
important dynamic processes. Maykut and Untersteiner (1969, 
1971) have constructed a detailed one-dimensional 
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thermodynamic model which provided a comprehensive view of the 
thermodynamic aspects. In the short term (days) , wind stress 
appears to be the dominant factor in ice prediction. On this 
time scale the thermodynamic effects of freezing and melting 
and the dynamic effects of weak ocean currents are small. As 
a result, predictions of sea-ice on this time scale can be 
reasonably accurate with only minimal consideration (simple 
parameterizations) of the ocean's influence. However, 
attempts to predict ice on a seasonal basis using either a 
thermodynamic or a dynamic model, or some combination of the 
two, without adequate consideration of the ocean's influence, 
have met with only limited success. The ice and ocean are 
closely linked over the long time scale and the interactions 
which occur between them must be accounted for. 

The U.S. Navy's current operational ice forecasting model 
run at the Fleet Numerical Oceanography Center (FNOC) is 
called the NORDA/FNOC Polar Ice Prediction System (PIPS) . One 
of the limitations of this model is its use of the monthly 
average ocean currents and ocean heat fluxes produced from the 
Hibler and Bryan (1987) ice-ocean model. This model was 
forced for several years of integration by repeatedly using 
a single year of observed atmospheric forcing data (December 
1978 to November 1979) . This year was chosen because it was 
the "FGGE" year in which a large number of drifting buoys were 
in the Arctic Basin which could provide a check for the model 
currents. The time period was not chosen for its similarity 
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to an average condition and, as noted by Hibler and Bryan 
(1987), it was in fact substantially different from the long 
term average. Additionally, the ocean below the mixed layer 
was constrained, with a three year damping time, to the 
observed circulation. Therefore, the ocean fluxes and 
currents calculated from the output of this model may very 
well be limiting the accuracy of the PIPS model predictions. 
This problem has yet to be resolved. 

The present work is aimed at improving our understanding 
and capability to numerically predict the location and 
thickness of Arctic ice on a seasonal time scale. This will 
be done through the further development of a linked ice-ocean 
numerical model using observed atmospheric forcing and a fully 
prognostic ocean model formulation. Several sensitivity 
studies will also be conducted to gain a better appreciation 
of the relative importance of the various mechanisms 
controlling the ice cover and to improve their representation 
in the model. The remainder of this chapter briefly describes 
previous efforts in ice-ocean modelling which have provided 
the background for this work. The specific goals and 
objectives of this research are also listed. Chapter II 
provides the geophysical context in which the model operates 
and Chapter III is a description of the model itself including 
the assumptions, boundary conditions and forcing fields under 
which it is run. Chapter IV examines the importance of 
including a fully interactive ocean model in studies of 
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interannual ice cover variations and Chapter V examines the 
effect of ice strength on those same variations and the ice 
thickness distribution. Chapter VI investigates the 
sensitivity of the model to frozen surface albedo changes and 
Chapter VII is a comprehensive examination of the various 
dynamic and thermodynamic processes included in the model and 
their relative importance to the evolution of the ice cover. 
Chapter VIII addresses the value of comparing the model 
simulations statistically and Chapter IX is a final summary 
which includes the conclusions from this work. 

B. PREVIOUS ARCTIC MODELLING 

Coupled ice-ocean models of the Arctic on seasonal time 
scales were developed only recently. Construction of such 
models was initiated primarily in response to the lack of 
quantified success of previous large scale ice modelling 
efforts by numerous investigators. General seasonal trends 
and basin-wide distributions similar to observed fields were 
achieved in these earlier efforts. However, the accuracy of 
regional ice distributions and/or the phase of the annual 
cycle were often poor. 

Manabe and Stouffer (1980) developed a coupled ice-ocean- 
atmosphere model which used a one-dimensional mixed layer 
ocean and a motionless slab with no leads for the ice. The 
atmospheric model provided daily variable forcing. This 
produced a wintertime ice cover about double that normally 
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observed in the Arctic. Washington et al . (1976, 1980) used 
a simple thermodynamic ice model but did not include any ice 
dynamics. Their model produced an excess of ice in the North 
Atlantic in all but the summer months and ice thicknesses in 
the central Arctic were less than half those interpreted from 
submarine surveys (LeSchack et al . , 1971; LeSchack, 1980). 

Once the emphasis shifted to prediction of the ice, the 
atmospheric portion of the coupled models was generally 
dropped in favor of prescribed atmospheric forcing and more 
elaborate ocean and ice models. The detailed thermodynamic 
ice model of Maykut and Untersteiner (1969, 1971) included the 
effects of ice salinity, brine pockets, shortwave radiation 
heating, vertical density variations, conductivity and 
specific heat of ice and water. However, this model required 
a great deal of computation to reach an equilibrium state. 
This made it impractical for application to a large scale 
three-dimensional coupled model. Semtner (1976a) simplified 
this complex model by parameterizing several of the processes. 
With the use of a three-layer version he was able to closely 
reproduce the results of Maykut and Untersteiner (1969,1971) 
in a wide variety of simulations, yet this simplified model 
had fewer computational requirements than one layer of ocean 
in a multi-level ocean model. This made it ideal for 
inclusion into a large three-dimensional gridded coupled 
model . 



8 



The Parkinson and Washington (1979) ice model used a 
further simplified "zero-layer" thermodynamic ice model (also 
proposed and tested successfully by Semtner (1976a)), in 
conjunction with a dynamic forcing model. The dynamics were 
based on stresses from wind, water, Coriolis force, internal 
ice resistance and ocean surface tilt. The atmosphere was 
represented by monthly climatological fields of temperature 
and pressure and the ocean was represented by a simple 30 m 
deep ocean mixed layer with a constant vertical heat flux. 
This model produced a reasonable yearly cycle of sea ice 
extent; however, the ice thickness distribution did not 
compare well with observational evidence. The combined 
thermodynamic-dynamic ice model concept was further improved 
by Hibler (1979, 1980). He also used the "zero-layer" 
thermodynamic approach of Semtner (1976a) but combined it with 
a more complex dynamic model and used daily observed 
atmospheric forcing data. The ocean was again represented by 
a simple fixed-depth mixed layer with a constant vertical heat 
flux applied beneath it. The essential difference in the 
Hibler approach was to link the dynamics of ice motion to the 
ice thickness by allowing the ice interaction to become 
stronger as the ice grew in thickness and/or concentration. 
In order to do this consistently, a viscous-plastic ice 
rheology was used. Arctic simulations with this model showed 
that the spatial distribution of ice thickness was strongly 
influenced by the ice dynamics. The results compared 
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favorably with the ice distribution observations of LeSchack 
(1980) and were a significant improvement over previous 
efforts. The improved ice dynamics as well as daily 
atmospheric forcing permitted a much better representation of 
lead formation. The proper representation of leads was 
considered important because heat and salt fluxes, as well as 
ice formation rates are nominally an order of magnitude larger 
in open leads than under the pack ice. 

Further analysis of the Semtner (1976a) thermodynamic ice 
model was conducted by Semtner (1984a). He concluded that an 
oversimplification of the thermodynamics tends to exaggerate 
the seasonal sea-ice cycle. The "three-layer" version of the 
thermodynamic model involves twice as much computation as the 
"zero-layer" version but does provide a much better 
representation of the ice. The computational requirements 
remain negligible compared to that of the ocean model. 
Therefore the "three-layer" model was recommended for use in 
large coupled models where accurate representation of the 
regional and temporal changes in the ice field are important. 

Several experiments have been conducted to determine the 
importance of atmospheric forcing. Three years of 
interannually variable forcing was applied by Hibler and Walsh 
(1982) to the Hibler (1979) ice model, but omitting inter- 
active ocean circulation effects. Their results again 
demonstrated the importance of ice dynamics for achieving 
accurate ice forecasts. Quantitative agreement with observed 
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interannually varying ice volume and transport was achieved; 
however, a noticeable seasonal bias was evident. The ice 
extent and thickness in winter and spring for all years was 
excessive while the extent and thickness of summer ice was too 
small with too much melting and excessive open water. They 
concluded that accurate parameterization of oceanic effects 
and a better treatment of the snow cover would probably be 
necessary in order for Arctic ice to be accurately modeled. 
Walsh et al. (1984, 1985) increased the interannually varying 
atmospheric forcing period to 30 years and used a more 
elaborate thermodynamic model with the Hibler (1979) ice 
dynamics. The seasonal bias in ice concentration was again 
evident in this model and it also showed considerable 
sensitivity to the number of vertical layers in the 
thermodynamic portion. 

An Arctic ocean circulation model was developed by 
Semtner (1976b) as an extension to the numerical ocean model 
of Bryan and Cox (Bryan, 1969) . Many of the observed features 
of the Arctic Ocean and Greenland Sea circulation were 
reproduced by this Arctic model despite the simplifications 
of mean annual wind forcing and a non-interactive ice cover. 
In order to develop an ice model suitable for seasonal ice 
predictions and daily variable atmospheric forcing, Hibler and 
Bryan (1984, 1986) coupled the Semtner (1976b) Arctic ocean 
model with a somewhat simplified Hibler (1979) thermodynamic- 
dynamic ice model. This resulted in the first seasonal, 
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coupled, ice-ocean numerical model. A linear term was 
included in the ocean model which damped the ocean's 
temperature and salinity values below the mixed layer to 
climatology. This ''robust-diagnostic” forcing was chosen to 
have a three year relaxation time. Multi-year drift from 
climatology was prevented; yet shorter term variations, 
particularly in the mixed layer, could still be prognostically 
simulated. The inclusion of the oceanic component greatly 
improved the simulated ice edge, removing much of the seasonal 
biases evident in the earlier work by Hibler and Walsh (1982) 
and Walsh et al . (1984, 1985) . Hibler and Bryan (1987) showed 
that many of the features of ice growth and decay, in various 
regions of the Arctic, result from a combination of numerous 
influences, often with opposing effects. For example, the 
East Greenland ice edge is forced southward by the East 
Greenland Current and by winter freezing. However, due to a 
strong vertical heat flux from the ocean, the amount of ice 
melted actually exceeds the amount frozen in this region. The 
southward extension of the ice zone is thus significantly 
inhibited . 

The Hibler and Bryan (1984, 1986) model simulated the 
seasonal variations in ice cover more accurately than any of 
the previous modelling efforts. However, Semtner (1987) noted 
that the three year diagnostic constraint may have dominated 
the forcing for the ocean model and thereby controlled the ice 
simulation. The constraint also prevented ocean circulation 
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and ice phenomena occurring on multi-year time scales (bottom 
water formation, carbon dioxide effects, etc.) from being 
forecast. Additionally, the climatological ocean was 
unrealistic in the variable-current, shallow-water regions of 
the Barents and East Siberian Seas. A fully prognostic, 
coupled, ice-ocean model would have no such limitations and 
was therefore developed (Semtner, 1987). This model, 
hereafter referred to as SM, is the basis of the current work. 
SM incorporated the best features of the previous modelling 
efforts. The Semtner (1976b) ocean model, optimized for high 
speed computers, was combined with the "three-layer" 
thermodynamic ice model (Semtner, 1976a) and a modified 
Hibler (1979) dynamic model. The daily forcing ice dynamics 
formulation of Hibler (1979) required more computation than 
a complete ocean model. This was not considered an equitable 
distribution of computation load. The Hibler dynamics were 
modified to use monthly averaged forcing instead of daily 
forcing, yet still retain the first-order dynamic effects 
(Hibler, 1988) . This also permitted a much longer time step 
and longer integrations. A simple 30 year mean seasonal cycle 
was computed from the forcing data used in Walsh et al. (1985) 
for atmospheric forcing. Monthly values of water inflow at 
the appropriate salinity, temperature and volume from rivers 
and the southern ocean boundaries were specified. These water 
mass fluxes were invariate in time. The model produced 
realistic representations of the average seasonal ice cycle 
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in distribution and concentration. The ocean circulation, 
including the major surface currents, water masses, and 
possible water transformation regions were also well simulated 
to the limits of the resolution. Results from SM emphasized 
the high degree of thermodynamic and dynamic interaction 
between the ice and ocean. 



C. OBJECTIVES 

The objectives of this work were to gain further 
understanding of the relative importance of mechanisms 
controlling the ice cover in the Arctic Basin and to develop 
a numerical model capable of accurately simulating the 
operationally and cl imatologically important interannual 
variations of ice cover in this region. In particular, the 
following points unique to this work were to be investigated: 

- Determine the importance of including a fully prognostic 
and interactive ocean component in a linked ice-ocean 
model. The importance was to be assessed relative to the 
accuracy of simulations of interannual variability in the 
Arctic ice cover. 

- Determine the sensitivity of the model to changes in the 
ice rheology strength parameter and subsequently select 
an optimum value of this parameter for incorporation into 
the model . 

- Determine the sensitivity of the model to changes of snow 
and ice albedo and if appropriate incorporate a more 
sophisticated albedo representation into the model. 

- Determine the relative importance of the various dynamic 
and thermodynamic ice cover forcing mechanisms included 
in the model, particularly the two features inherent in 
the ocean model (vertical ocean heat flux and ocean 
current stress) . 
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- Assess the regional dependency of the above ice cover 
control mechanisms. 

- Assess the statistical significance of the differences 
between the various model run results. 

D. METHODOLOGY 

This research was accomplished through the utilization and 
further development of the coupled numerical model of Semtner 
(1987). The Semtner (1987) model is the most freely 
prognostic, physically representative and computationally 
efficient of the large-scale, linked, ice-ocean numerical 
models previously developed for the Arctic. Considerable 
skill has been demonstrated with the model in reproducing the 
annual cycle of ice extent with climatological forcing. Most 
importantly, the inclusion of an interactive ocean reduced the 
pronounced tendency of previous ice-only simulations to 
exaggerate the annual ice cycle. 

Numerous model integrations were conducted. In each 
experiment the model was modified to examine some particular 
mechanism which was expected to influence the interannual 
variability of the ice cover. In the majority of cases the 
model was integrated for ten years. However, some of the 
sensitivity testing required a large number of repeated runs 
and in those cases the model runs were limited to one 
integration year. Additionally, in the sensitivity runs where 
inclusion of an interactive ocean was of limited concern, the 
model was run with a prescribed ten-year average ocean 
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condition vice a prognostic ocean component. The ice cover 
simulations were evaluated from January 1971 to December 1980 
inclusive. These years were chosen because the observed ice 
concentration data set, to which the modelled ice 
concentrations were compared, was most accurate in this 
period. This was especially true after 1972 due to the 
assimilation of microwave remotely sensed data. Parkinson et 
al. (1987) was also an excellent source of comparison data for 
the years 1973-1976. 

The observed ice concentration fields were provided by 
Walsh, and were the same as those used in Fleming (1987) but 
converted through bi-linear interpolation to the grid used 
here. The observed ice concentration fields were regarded as 
the "true" monthly ice concentrations. However, they were in 
fact the average of data received from several sources, each 
with its own source of error (Walsh and Johnson, 1979) . 
Nevertheless, these fields were still the most comprehensive 
and accurate representations of real data available and 
covered a long enough time period to permit viable statistical 
analysis. All 120 months of observed ice concentration were 
contoured to gain a qualitative appreciation of the annual 
cycle, its variability and absolute extent. This also 
provided the baseline view for visual comparison with the 
simulated ice concentration data. 

The seasonal average ice draft contours from Bourke and 
Garrett (1987) were the primary source of comparison data for 
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the simulated ice thicknesses and thickness distributions. 
Although these contours represent "observed” data from 
submarine transits, the data were heavily averaged in time and 
space to produce the contours. Nevertheless, the general 
thickness and distribution patterns did provide an indication 
of the true seasonal ice thickness distribution. 

The monthly varying atmospheric forcing fields were also 
provided by Walsh. These fields were produced from the same 
data used by Walsh et al. (1985) and were also the source of 
the 30-year mean forcing using in Semtner (1987). 

It was shown by Fleming (1987) and Garcia (1988) that 
correlations between ice concentration and thermodynamic 
forcing variables (sea surface temperature (SST) and air 
temperature) and dynamic forcing variables (winds and previous 
ice condition) vary significantly between regions which 
experience different dynamic ocean conditions. To gain some 
appreciation for these regional differences but still maintain 
the amount of data manipulation within reasonable limits, the 
complete grid space was divided into four regions (see Figure 
1.1). The divisions were selected based on the general ocean 
current regime. The four regimes represented were: 

- Region 1 X=l,18 Beaufort Gyre. 

- Region 2 X=19,27 Transpolar Drift, Ellesmere Is. 

and Greenland ice convergence 
zone . 
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Figure 1.1 Geography of the Arctic Ocean as represented in 
the model. The small boxes indicate the centers 
of land gridpoints. Latitude circles every two 
degrees are shown over the ocean gridpoints. The 
arrows point to the mouths of the indicated rivers 
and the numbers indicate the regions (from 
Semtner, 1987) . 
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- Region 3 X=28,34 Transpolar Drift (TPD) , East 

Greeenland Current (EGC) . 

- Region 4 X=35,45 Norwegian Current. 

where X defines the location of the region on the X axis of 
the grid. Each region extends from 1 to 42 on the Y axis of 
the grid. In those cases where the modified model simulations 
produced ice fields more representative of the observed ice 
fields than before the modifications were made, the 
appropriate changes were incorporated into the model. 
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II. ARCTIC CLIMATOLOGY 



This chapter presents the geophysical characteristics of 
the Arctic region emphasizing those features which are most 
relevant to numerical modelling of the ice cover and ocean. 

The primary dynamic driving forces for sea ice growth and 
movement are the air stress resulting from wind-induced 
surface drag and the water stress resulting from ocean 
currents. Variations in wind forcing are determined by the 
passage of synoptic scale weather systems having spatial 
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fluctuations with length scales of approximately 10 km and 
time scales of hours to days. Inertial oscillations of 
drifting sea ice have similar periods. Longer term 
fluctuations, with time scales of weeks to seasons, represent 
responses to atmospheric and ocean forcing integrated over 
equivalent time scales. No firm evidence exists to suggest 
that small time and spatial forcing is unimportant to the 
seasonal evolution of the ice cover. However numerical 
modelling of the Arctic ice has not developed sufficiently at 
present to resolve such forcing for long periods over the 
entire Arctic basin. Therefore, the emphasis in this work 
will be on the monthly time and larger spatial scale features. 
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When examining fluctuations varying on time scales of a 
season to several decades, the dominant signal in sea ice 
variability is clearly the annual cycle (Parkinson et al., 
1987). Figure 2.1 shows the 25-year (1953-1977) envelopes of 
the positions of the Arctic ice edge at the end of August and 
February. At most longitudes the seasonal change from the 
summer to the winter ice edge positions is considerably 
greater than the 25-year range of extremes for a particular 
month. The overall annual freeze-melt cycle is primarily a 
response to the seasonally varying insolation, albeit with a 
large amount of feedback and interaction with the ocean. 



A. ATMOSPHERE 

The climate in the Arctic region is largely determined by 
four aspects (Sater et al., 1971): 

- A distinctive cycle of prolonged periods of daylight with 
net positive surface heat balance followed by prolonged 
periods of darkness with a net negative surface heat 
balance . 

- A very persistent, cold-cored, circumpolar, upper level 
vortex which steers the mesoscale surface-level weather 
systems . 

- A surface cover of snow or ice for periods of the year, 
which causes significant changes in surface albedo and 
therefore surface heat absorption. 

- A strong temperature inversion above snow or ice surfaces 
which acts to decouple the surface from higher level 
atmospheric activity. 

Many factors are involved in determining the net surface 
heat balance for any particular region of the Arctic. These 
will not all be examined here. Suffice it to say that the 
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Figure 2 . 1 




Composite maximum and minimum ice extents at the 
end of August and February. The hatched region 
is the difference between extreme positions of the 
0.5 ice concentration lines based on 25 August 
grids and 25 February grids (1953-1977) . The dots 
represent the grid points used in digitization of 
monthly ice concentration data (from Walsh and 
Johnson, 1979) . 
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main factor is believed to be the insolation rate. Figure 2.2 
summarizes the effects in a radiation balance diagram 
determined for the Barrow Strait region. The diagram shows 
long periods of net heat loss which correspond to the polar 
winter and periods of large ice production. The long heat 
loss periods are followed by short intense periods of net heat 
gain corresponding to the long daylight hours and the ice 
melting of summer. This figure is reasonably representative 
of the heat balance in the Arctic in general. However, a term 
which is missing and which is regionally and spatially quite 
variable is the advective heat flux, Q v . This term represents 
the heat supplied to the surface from vertical and horizontal 
advection in the ocean. Mixing processes could also be 
included. Huyer and Barber (1970) suggested that this term 
was negligible in the Barrow Strait. However, in other 
regions of the Arctic, large differences in latitudinal extent 
of the ice exist as can be seen in Figure 2.1. This suggests 
that Q v may be very significant to the surface heat balance in 
other regions in the Arctic. 

Surface weather systems, especially cyclones, are 
primarily generated in regions of strong differential heating. 
Significant differences in surface temperatures occur at land- 
sea, ice-sea, and warm-cold current boundaries. Several such 
boundaries are found in the Arctic in fairly constant 
positions. Weather systems are therefore more prevalent in 
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Figure 2.2 The radiation balance showing the short-wave 
incident radiation, Q, the absorbed short-wave 
radiation, Q a , the effective back radiation, Q b , 
and the net radiation, Q n (from Huyer and Barber, 
1970) . 



these regions, but their subsequent influence on other areas 
is dependent on the upper level steering and the conditions 
of the surface over which they pass. The cyclones will tend 
to increase in intensity if they remain over relatively warm 
open water, which provides heat and moisture. However, they 
will tend to weaken and occlude if they travel over cold and 
dry (ice or land) surfaces. These factors imply that 
considerable variability exists in the dynamic atmospheric 
forcing between different regions and in different seasons in 
the Arctic. Climatological averages or even monthly averages 
of atmospheric forcing will undoubtedly mask some of this 
variability. 

Variations in surface albedo cause the percentage of solar 
radiation absorbed to differ by as much as 7 0% between 
different types of snow, ice and open water. Robinson et al., 
(1986) have used comprehensive sets of satellite-derived data 
and ground confirmation data to identify four ice surface 
albedo classes in the Arctic. Albedos range from 0.80 for a 
surface with fresh snow covering 95% of the ice to 0.29 for 
a heavily ponded surface with less that 10% bare snow or ice. 
Kukla and Robinson (1980) determined that 0.10 was a 
reasonable albedo for open water in the high latitudes. These 
differences cause the melting and freezing processes to 
accelerate rapidly once started. Changes in the ice cover can 
be very dramatic over short periods of time. The Arctic 
stratus deck, common to the region in summer, also influences 
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surface albedo significantly, but this effect is very 
difficult to quantify. Considerable variability exists in the 
literature as to the correct albedo to use in large scale 
modeling. Numerous authors have shown that numerical 
thermodynamic ice models can be very sensitive to this 
variable (e.g., Shine and Henderson-Sel lers , 1985; Ross and 
Walsh, 1987; Semtner, 1976a). Inclusion of a realistic albedo 
representation is therefore difficult but probably important 
in an accurate numerical sea-ice model. 

Open leads probably account for less than 10% of the 
surface area encompassed by the Arctic ice; however the heat 
and moisture exchanges which occur above them can be an order 
of magnitude higher due to the large differences in 
temperature and radiation characteristics between seawater and 
the pack ice surface. This results in localized areas of high 
static instability, cloud formation and large heat and 
moisture fluxes. Leads are generally too small to resolve 
using current large-scale numerical ice models; however, they 
are important enough to require that their effects be 
parameterized. 

The strong temperature inversion common to much of the 
Arctic is the result of the very cold ice and land surfaces. 
In winter the inversion layer can deepen to the 850 mb level, 
becoming most intense under calm, clear anti-cyclonic 
conditions (Sater et al . , 1971). The inversion is a very 
persistent feature which will quickly reform after periods of 
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strong winds, cloud cover or storm precipitation as long as 
surface temperatures remain relatively low. In summer this 
feature is confined to the polar ice cap as surface 
temperatures over the ice-free water at more southern 
latitudes become too high. The inversion acts- to insulate or 
decouple the surface from upper level winds and also tends to 
put a cap on any moisture which may be absorbed in the 
boundary layer. The result is relatively calm (average less 
than 5 m/s) surface wind conditions and a predominance of 
stratus clouds, especially when the melt season commences. 
Care must be taken to ensure that surface winds calculated 
from higher level pressure charts, if used as forcing fields, 
take this feature into account to ensure realistic wind 
forcing on the ice. 

B. WATER MASSES 

The waters of the Arctic Seas are often described on the 
basis of temperature and salinity. As such, they are 
comprised of three main water masses: Arctic Surface Water, 
Intermediate or Atlantic Water, and Deep or Bottom Water 
(Coachman and Aagaard, 1974). Arctic surface water is 
generally limited in depth to about 200 m. It has the most 
variable characteristics and can be modified by the weather, 
the season, and/or the physical environment. Temperatures in 
this layer vary from -1°C to over 2°C. The salinity may be 
uniform to approximately 50 m, below which a sharp halocline 



27 



increases the salinity to about 34.5 ppt at the bottom of the 
layer. The variety of conditions in the surface water is 
evident in Figure 2.3, where vertical profiles of temperature 
and salinity from a variety of Arctic basins are plotted. 
Coachman and Aagaard (1974, p. 9) state that: 

The most important processes conditioning and 
modifying the surface layer are: 

- Addition of mass (fresh water) from the land, 
primarily from the large Siberian rivers; 

- Addition of fresh water locally through melting of 
ice ; 

- Heat gain through absorption of solar radiation in 
non-ice-covered areas during summer; 

- Concentration of salt and hence increase of density 
of surface water, through freezing of ice; 

- Heat loss to the atmosphere through any open water 
surface, including leads in the central Arctic pack 
ice ; and 

- Inflow and subsequent mixing of Atlantic and Pacific 
waters . 

Processes 1, 2, and 3 normally occur only from June to 

September and lead to a decrease in water density. These 
buoyant waters form a surface cap which absorbs radiation and 
warms. Therefore, in summer, ice free regions tend to have 
warmer and less saline surface layers. In areas where the ice 
does not recede, surface temperatures remain near freezing as 
incoming energy is used to melt the ice, but surface layer 
salinities are reduced due to ice melting. 
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Figure 2.3 



Vertical profiles of temperature and salinity for 
various northern high-latitude basins (from 
Coachman and Aagaard, 1974) . 
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Processes 4 and 5 have the greatest impact in winter. In 
some areas, such as the shelf waters, conditions are such that 
the water becomes dense enough to penetrate into the 
intermediate Atlantic layer or in extreme cases even form 
bottom water (Aagaard et al., 1985) . However, in general, the 
strong pycnocline at the base of the surface layer prevents 
mixing due to surface density changes from penetrating below 
200 m (Coachman and Aagaard, 1974) . 

As a consequence of the halocline-derived mixing barrier, 
process 6 does not exert a great influence on the surface 
water variability. However, regions of high surface 
salinities (33 to 34.5 ppt) are found in the Greenland and 
Labrador Seas and east Baffin Bay. These high salinities are 
due to the advection of North Atlantic surface water into the 
Arctic via warm surface currents from the south. One further 
consequence of the halocline barrier is that the ice cover is 
often insulated from the relatively warm Atlantic layer. 
However, in some locations like the continental slopes, the 
Atlantic layer is forced to shallower-than-normal depths. 
Vigorous surface mixing can then break through the halocline 
and vertical heat fluxes can occur (Coachman and Aagaard, 
1974) . This results in the ice melting from the bottom or not 
forming at all, increasing the area of open water. 

Many mixing events are small in scale but complicated and 
intense. In a large scale ice model these effects should be 
included, but as was the case for leads, they need to be 
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parameterized in some manner due to the scale limitations of 
the model . 

The second water mass is called the Atlantic layer. This 
warm and salty (35.0 to 35.1 ppt) water originates as North 
Atlantic water flowing into the Arctic basin through the Fram 
Strait and the Barents Sea (Weigel, 1987). It extends from 
200 m to 900 m with temperatures above 0°C. A temperature 
maximum of approximately 0.45°C, observed throughout the 
central and western Arctic Basin, occurs between 300 m and 500 
m depth, dependent on location. The salinity gradually 
decreases to approximately 34.9 ppt in the Arctic Ocean and 
Greenland Sea and approximately 34.6 ppt in Baffin Bay. The 
intermediate layer is a source of heat which can be tapped by 
vertical circulation. As a result it may play an important 
role in the net heat balance of the surface mixed layer. The 
net heat balance in the mixed layer ultimately determines 
whether or not ice will form. The ocean model used here has 
been shown to simulate this warm intermediate layer quite well 
(Semtner, 1984, see Figure 2.4). 

The third water mass, the deep and bottom waters, have 
temperatures below 0°C. The salinity is nearly constant from 
the bottom of the Atlantic layer to the ocean floor. The 
intermediate Atlantic water and deep water are advected into 
and out of the Arctic seas from adjacent areas, principally 
through Fram Strait in the North Atlantic. Both water masses 
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Figure 2.4 Vertical sections of observed temperature and 
salinity along the prime meridian, from Coachman 
and Aagaard (1974) (upper panels) ; and vertical 
sections of simulated temperature and salinity 
across the model domain at approximately 20 
degrees E (lower panels). The lower panels have 
a stretched vertical scale, (from Semtner, 1984b) 
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are nearly isohaline and therefore of nearly uniform potential 
density. 

A temperature difference exists between the deep and 
bottom waters of the Canadian and Eurasian Basins. The deep 
Canadian Basin averages approximately -0.45°C while the 
Eurasian basin temperature is approximately -0.79°C (Aagaard 
et al., 1985). The deep and bottom waters of the two basins 
are kept separated by the Lomonosov Ridge which acts as a sill 
between them. 

Observations of the vertical structure at various stations 
have been taken over long time periods and in many different 
locations. A remarkable similarity in the profiles has led 
to the conclusion that the Arctic basins are in a long-term 
dynamic steady-state condition (Coachman and Barnes, 1961) . 
It was further noted that observed distributions of Arctic 
water properties are a result of continuing processes within 
the basins. Therefore, surface water temperature-salinity (T- 
S) profiles reflect the local modifying processes, while T-S 
profiles for depths below 200 m reflect the common origin of 
the water. 

The net heat balance ultimately determines how much ice 
is produced or melted. The heat flux provided by the ocean 
is an important parameter in determining the net heat balance 
at each grid point as noted in Chapter I. The ocean component 
of a linked ice-ocean model should therefore have sufficient 
resolution to reasonably define the boundaries of the various 
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water masses in three dimensions. It must also be able to 



simulate the water motion which advects oceanic heat to the 
ice and convects the cold saline water to the bottom. 

C. CURRENTS 

The surface circulation of the Arctic Basin has been 
derived from satellite observations, the trajectory of ice 
islands, buoys and floe stations, and the analysis of ship's 
tracks. The circulation of the Arctic waters is due both to 
water density differences and wind forcing. Large anomalies 
in the flow (compared to the long-term mean currents) are 
often observed; however, the long-term mean surface currents 
have been reasonably well documented in the northern Atlantic 
and Arctic Seas (Krauss, 1986). Figure 2.5 covers the area 
of interest in this study and indicates the major ocean 
currents. Recent moored current measurements have allowed a 
reassessment of the boundary and interior circulation in the 
Arctic Basin (Aagaard, 1988) . Aagaard emphasizes the 
important role of mesoscale eddies in the Canadian Basin and 
proposes a separate cyclonic flow for the Eurasian Basin. 
This newest assessment is very similar to the simulated 
Atlantic layer circulation in Semtner (1987) and is a 
significant change from the previous basin-wide cyclonic 
circulation proposed by Coachman and Aagaard (1974) . Figure 
2.6 shows this comparison. 



34 




Figure 2.5 Schematic of the large-scale horizontal 
circulation patterns in the surface waters of the 
Arctic region (from Parkinson et al., 1987). 
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Figure 2.6 Comparison of the simulated circulation in the 
Atlantic layer (565 m) (upper diagram) (from 
Semtner, 1987) with the reassessed subsurface 
circulation presented by Aagaard (lower diagram) 
(from Aagaard, 1988). 
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The Arctic Ocean contains two main surface flow patterns. 
In the Canadian Basin, the Beaufort Sea gyre rotates 
clockwise. It is driven primarily by the mean wind. Across 
the Lomonosov Ridge, the Transpolar Drift Stream (TPD) flows 
directly across the Eurasian Basin, from the Laptev Shelf to 
the Fram Strait where it exits as the East Greenland Current 
(EGC) . Several small, weak gyres split off the Transpolar 
Drift Stream as it passes close to the various islands off the 
coast of the USSR. 

Circulation and ice movement patterns in the peripheral 
Atlantic seas of the Arctic Ocean are dominated by two major 
currents systems, one warm, the other cold (Krauss, 1986) . 
The cold currents support ice growth and enhance the spatial 
extent of sea ice, while the warm currents advect heat and 
melt the ice or preclude its formation. The southerly extent 
of the ice edge between regions of dissimilar temperatures can 
vary as much as 30 degrees of latitude. 

Sea ice off the east coast of Greenland is found to 
originate primarily in the Arctic basin and as a result can 
be quite thick. The EGC continually carries a wide belt of 
Arctic pack ice southward through Fram Strait. As the ice 
season develops, the ice edge extends further south down the 
east Greenland coast. In extreme years the ice edge will 
extend as far south as Kap Farvel or as far north as 7 0 
degrees N (Sater et al., 1971). 
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Strong contrasts in ice coverage are also observed within 
the Baffin Bay-Davis Strait region. The warm West Greenland 
Current (WGC) follows the bathymetric contours of the southern 
and western Greenland continental shelf, keeping the southwest 
coast of Greenland ice free during most winters. However, the 
cold Labrador Current contributes to the large concentrations 
of ice found along Canada's east coast, well south of Baffin 
Island. This current also carries the majority of icebergs 
southward through the Grand Banks fishing zones, the Hibernia 
oil fields and into the North Atlantic sea lanes (CIA, 1978) . 

A reasonably accurate representation of the ocean currents 
is obviously important if sea ice is to be accurately 
simulated. The currents provide a direct influence through 
dynamic forcing as well as a thermodynamic influence via heat 
advection. 

D. ICE 

Sea ice is a feature which contributes to the uniqueness 
of the oceanography in the Arctic region. Coachman and 

Aagaard (1974) state: 

The general oceanographic consequences of a perennial 
or seasonal ice cover are: 

- The water temperature of the near-surface layer in the 
presence of ice is always maintained close to the 
freezing point for its salinity by the change of phase 
process . 

- Salt is excluded from the ice to a varying extent, 
but the water under the ice is always enriched in salt 
by any ice growth. The dependence of water density 
on temperature and salinity is such that close to the 
freezing point, density is almost solely a function 
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of salinity. Therefore, ice formation can increase 
the density locally and some vertical convection may 
result . 

- In the transfer of momentum from the atmosphere to the 
ocean, the wind must act on the sea through the 
intermediary of the ice. 

The Arctic Ocean ice pack is confined by a nearly 
continuous boundary of land. The associated constraint on 
equatorward transport is a major reason why Arctic ice 
survives longer and develops into more complex forms than ice 
found in southern polar regions. An annual net heat loss and 
stratification of the underlying water also contribute to the 
longevity of Arctic ice. 

Ice formation in open water starts in autumn. As days 
grow shorter and nights longer, the amount of solar insolation 
decreases. Since long-wave radiation from the Arctic remains 
approximately constant during the year, the energy budget 
changes from a net gain in summer through equilibrium to a net 
loss in the fall (Figure 2.2). The relatively warm mixed 
ocean surface layer of summer is cooled until it reaches its 
freezing point and ice crystal formation commences. Ice 
initially forms around the boundaries of the polar ice pack 
and over the shallow protected waters of high-latitude 
coastlines. The ice-covered areas continue to expand until 
they merge, forming an ice-locked Arctic Ocean from October 
to June. The ice cover starts from a summer minimum of 
approximately 5.2 x 10 6 km 2 then more than doubles in areal 
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extent to a maximum of 11.7 x 10 6 km 2 by the end of the ice 
season (CIA, 1978, p. 12). 

The ice cover grows continuously until spring. Once the 
days start to lengthen, the amount of solar radiation 
increases until the surface energy budget is again positive. 
As noted earlier, snow cover will reflect approximately 80 
percent of the sun's incident radiation, which slows the 
initial heating stages. However, once the air temperature 
reaches the snow's melting point, the albedo rapidly drops to 
50 percent and lower. This results in a period of rapid snow 
melt. Continued heating initiates melting of the underlying 
ice, causing cracks and flaws to develop. The surface melt 
water drains through the cracks, further eroding them until 
the ice breaks up into floes. Eventually, the ocean surface 
layer heats up, the floes gradually dissolve, and the cycle 
is complete (Langleben, 1972) . 

Extremes of ice cover in summer and winter for the period 
1953-1977 were shown in Figure 2.1. Obviously, the sea ice 
cover experiences a large degree of variability both 
seasonally and interannually . Ice cover is something of a 
misnomer as the sea surface is rarely covered by an unbroken 
expanse of ice. Even in the very thick ice regions of the 
winter polar pack, infrared measurements indicate that up to 
ten percent of the area of the Arctic ocean consists of open 
water or thin ice from recently refrozen leads (Sater et al., 
1971, p. 41) . 
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Multi-year ice (ice which has survived a summer's melting) 
comprises the majority of the polar pack, which averages 3 m 
in depth (Bourke and Garret, 1987) . First-year ice rarely 
grows to a thickness greater than 2 m. However, the depth of 
the ice in any location is largely dependent on the external 
forces of wind and currents. These cause the ice to converge 
and diverge. When a region of ice converges, it buckles, 
folds, and overlaps, forming a rugged terrain and areas of 
considerably thick ice. For example, the Beaufort Gyre's 
anticyclonic flow causes ice to converge along the north coast 
of Ellesmere Island and Greenland. The number of ridges in 
this region is well above the average for the Arctic pack 
(Weeks et al., 1971). The mean ice thickness here is of the 
order of 6 m to 8 m (Bourke and Garrett, 1987). In contrast, 
the ice pack east of Spitsbergen (Svalbard) is not confined 
by land and is free to diverge. Here, the average ice 
thickness is significantly less, averaging approximately 2 m. 

The thickness, distribution and concentration of sea-ice 
are influenced by both dynamic and thermodynamic factors. The 
ice portion of an accurate ice-ocean linked model needs to 
include as much of the physics of these mechanisms as 
possible. 
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COUPLED NUMERICAL MODEL 



The coupled model used in this study contains an ice model 
with both dynamic and thermodynamic forcing linked to a fully 
prognostic primitive equation ocean model. The grid domain 
includes the Norwegian and Greenland Seas, and the entire 
Arctic basin. The geography is smoothed due to the 
limitations of the resolution leaving only Spitsbergen and 
Iceland as true islands. Bottom topography is included. The 
model is forced by wind stress and energy balances at the sea 
surface and water mass balance at open boundaries. The ocean 
model is coupled with the ice model by momentum, heat, and 
salt exchanges through a common ocean mixed layer. Much of 
this chapter follows the very clear descriptions of Van 
Ypersele (1986) who applied a similar model to the Antarctic. 

A. OCEAN MODEL 

The physical state of the ocean is characterized by seven 
quantities: pressure, potential temperature, salinity, 
density, and the three components of velocity. A three- 
dimensional ocean model must be able to predict the evolution 
of those seven quantities from an initial state when adequate 
boundary conditions are specified. 

Seven equations are required to relate the seven variables 
and close the system. These equations are derived from basic 
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principles including the conservation of momentum, energy and 
mass as well as Newton's law of gravitation. An equation of 
state is used to calculate density from temperature, salinity, 
and depth. The equations must be simplified yet still retain 
the most important terms in order to solve them numerically 
and produce reasonable results. Following Semtner (1974, 
1986b) the necessary simplification was provided by making the 
following approximations and assumptions about the system of 
equations: 

- The ocean is incompressible. This eliminates acoustic 
waves from the possible solutions and leads to a simpler 
continuity equation without density terms. It is a 
reasonable assumption because the volume changes due to 
pressure are small. For example, the change of pressure 
corresponding to a change of depth of 1000 m would alter 
the volume of a typical sample of seawater by less than 
0.5%. (Pond and Pickard, 1983). 

- The ocean behaves as a horizontally isotropic fluid with 
constant kinematic eddy-viscosity coefficients in the 
vertical and horizontal directions. This means that an 
error is introduced in the mixing of momentum if the 
isopycnals are not horizontal. 

- Exchanges of heat and salt which occur at sub-grid scale 
can be represented by similar eddy-dif fusivity 
coefficients in the vertical and horizontal directions. 
On the molecular scale, salt diffusion in water is on the 
order of 100 times slower than heat diffusion. However, 
turbulent processes tend to dominate which permits use of 
the same eddy-dif fusivity coefficients for both salt and 
heat . 

- The thin-shell approximation is made. This allows certain 
terms in the horizontal momentum equations to be neglected 
because the depth of the ocean is negligible compared to 
the radius of the earth. 

- The hydrostatic assumption is made. The differential of 
pressure with respect to depth is considered only a 
function of the product of density and the acceleration 
of gravity. This approximation is a result of scale 
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analysis of the vertical momentum equation which reduces 
to the form: 




= ~Pg 



(3.1) 



- The Boussinesq approximation is made. This is basically 
the neglect of small density variation effects on the 
fluid momentum and therefore on the acceleration for a 
given force, while maintaining these same effects on 
weight in the earth's gravity field (buoyancy) . 
Consequently, average density over the domain can be used 
in the horizontal momentum equations. However, the actual 
density values are required to derive pressure from 
Equation 3.1. 

- The Coriolis terms and the viscous terms involving the 
vertical speed, w, in the horizontal momentum equations 
are neglected on the basis of scale analysis. 

- The rigid lid approximation is made. This is achieved by 
specifying that the vertical velocity, w, be zero at the 
upper surface. This process filters out high-speed 
surface gravity waves in order to permit use of a longer 
timestep . 

The three momentum equations are derived from Newton's 
second law of motion. The resultant equations in vector form 
are often called the Navier-Stokes equations. They provide 
a relation between the velocity, the pressure gradient force, 
the viscous stresses, and the body forces such as gravity and 
Coriolis (e.g., Tritton, 1977). 

After the simplifications are made, as allowed by the 
assumptions and approximations described above, the momentum 
equations in spherical coordinates become: 
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where f is the Coriolis parameter (f = 2flsin0 [ s’ 1 ] ) , fi is the 
angular speed of rotation of the earth, (fl = 7.272 x 10' 5 
rad/s) . A z and A h are the vertical and horizontal eddy 
viscosity coefficients, respectively (0.3 x 10' 4 m 2 /s and 12 
x 10 4 m 2 /s) . r e is the earth's radius (6.37 x 10 8 m) . L(cr) 
is the advection operator which in spherical coordinates 
becomes : 



L(a) ~ 7T^0 JfX ( ua ) + F5o^o % (vacos ^ + & (w«) 

e e 
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The horizontal Laplacian operator becomes: 
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e e 

where a is an arbitrary quantity. 

Equation 3.4 is the hydrostatic equation governing the 
depth related variation of pressure in the ocean. It assumes 
that a perfect balance of forces exists in the vertical (i.e., 






(3.6) 
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vertical accelerations are negligible) . This is valid for 
large scale motions, but it is not adequate to describe small- 
scale convection, where horizontal and vertical scales of 



motion tend to be similar. Therefore, the model must be 
supplemented by an implicit treatment of convection, which 
will prevent the unrealistic condition of unstable 
stratification. The convective adjustment scheme adopted here 
is such that if the condition 



is verified, then the new values of temperature and salinity 
are recomputed from 



Level k + 1 is taken to be the common reference level, and ( — ) 
indicates a weighted average over layers k and k + 1. 

The conservation of mass, together with the assumption of 
incompressibility, gives the continuity equation: 



In spherical coordinates this becomes: 
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The theory of thermal transfer in a fluid of constant heat 
capacity, together with an eddy-dif fusivity closure 
hypothesis, gives a conservation equation relating the 
evolution of temperature to the motion field and the first 
and second spatial derivatives of the potential temperature: 



§ + L(0) = + KhVf,0 (3.1D 

K z and K h are the vertical and horizontal eddy diffusivity 
coefficients, respectively (0.3 x 10' 4 m 2 /s and 2 x 10 3 m 2 /s) . 
A similar equation is derived for salinity: 
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An equation of state is used to express the density as a 
function of potential temperature, salinity, and pressure: 



P = p(B,S,p) 



(3.13) 



For the purposes of numerical modelling, the equation of 
state is expressed as a polynomial approximation to the 
results of laboratory experiments. The approximation of 
Eckart (1958), utilized by Semtner (1974) will also be used 
for this study. 
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B. ICE MODEL 



The sea-ice portion of the coupled model contains both 
thermodynamics and dynamics. The thermodynamic portion is 
based on the three-layer thermodynamic ice model proposed in 
Semtner (1976a) , which is a simplification of the 
comprehensive ice model of Maykut and Untersteiner (1971) . 
It is forced from above by the atmospheric fluxes computed 
from atmospheric data, and from below by oceanic heat fluxes 
prescribed by the ocean model. The dynamic portion is a 
simplification of the "viscous-plastic" rheology approach used 
by Hibler (1979). The simplification reduces the 
computational load significantly and permits use of the same 
monthly averaged atmospheric data for forcing as used by the 
ocean model. 

For the three-layer thermodynamic model, sea-ice is 
assumed to be a horizontally uniform slab of ice, represented 
by two layers of equal thickness. A covering layer of snow 
is possible. When ice is present, the temperature in the 
oceanic mixed layer below the ice cover is assumed to be at 
the freezing point of seawater (-1.9°C). The temperature 
profile through the ice and snow is governed by the one- 
dimensional heat equation: 



r UY _ . &T 
pL, W ~ k c^ 



(3.14) 
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where p is the density [900 and 330 kg in" 3 for ice and snow, 
respectively], C is the specific heat capacity [2.09 kJ kg' 1 
K' 1 for both ice and snow], T is the temperature [K], and k 
is the thermal conductivity [4.068 and 0.31 W m^K" 1 for ice 
and snow, respectively] . 

At the top and bottom surfaces of the cover (ice or snow 
as applicable) , a balance of fluxes is forced to exist such 
that any heat excess or deficit that occurs at either boundary 
is used to melt or freeze, respectively. Melting can occur 
at either surface; however freezing can only occur at the ice 
bottom. Therefore when a negative flux imbalance occurs at 
the surface, the heat equation is applied to the ice to 
ascertain the effect at the ice bottom before the amount of 
freezing is determined. The melting point for snow is set at 
0°C and for ice -0.1°C. The heat of fusion for snow is set at 
110 MJ/m 3 and for ice at 301 MJ/m 3 . Diffusive fluxes must be 
equal at the internal snow-ice interface. 

By late summer, the ice will have completely melted in 
certain regions. A heat budget equation is then applied to 
the vertically isothermal surface mixed layer and its 
temperature is allowed to increase. Horizontal mixing and 
advection of the warmed oceanic surface layer can then spread 
the heat under adjacent ice covered areas to speed the melting 
process there. Once the freezing season commences, the heat 
budget cools the mixed layer to its freezing point, ice starts 
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to form and the heat equation for ice/snow must again be 
applied. 

In order to determine the final temperature of the upper 
surface when ice is present, several conditions are imposed 
to control the way the various fluxes are applied. If the 
initial temperature at the upper surface, T surf < T,,^ , the 

upper surface flux balance equation is given by: 

Q S w( 1_I °) (!— a surf) + Q[w “ Ql\V + ^Ii + ^LE + ^cond = 0 

(3.15a) 

where a surf is the surface albedo and I 0 is the fraction of 
solar radiation absorbed in the ice cover. I 0 is set equal 
to 0 for snow and to 0.17 for bare ice. The remainder of the 
terms are: 

- Q lu longwave heat flux, the arrow indicates downward (+ve) 
or upward (-ve) . 

- Q h sensible heat flux (+ve downward) . 

- Q le latent heat flux. 

- Q su short wave solar flux. 

- Q cond flux conducted downward from the surface to the 
bottom of the ice. 

These are described in more detail in section E (Forcing 
Fields) . 

If T $urf = T me(t , the upper surface flux balance equation is 
given by: 
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(3 . 15b) 



Qsw( 1-I °) ^~ asurf ) + Qlw “ ^LW + % + ^LE + ^ C0lld 



-L 




z~0 



where hj and h s are the thicknesses of the ice and snow, 
respectively, and Ly, is the latent heat of fusion. 

At the bottom of the ice, the flux balance equation is: 



Q 



cond 



-Q -I n - dh I 

ocean - 



z=h; 



(3 . 16) 



where Q ocean is the oceanic heat flux [W m' 2 ] . 

The outgoing longwave radiation Q LU follows the Stefan- 
Boltzmann law: 



Q 




to T-1 

surf 



(3.17) 



where e is the surface emissivity (non-dimensional, and 
assumed equal to 1) , and a is the Stef an-Boltzmann constant 
[= 5.67 x 10' 8 W m' 2 K' 4 ] . 

The final temperature at the upper surface, T surf , is 
obtained by linearizing the blackbody emission term 3.17. 

The conductive flux Q cond inside the snow and ice is 
computed by assuming a linear vertical temperature profile 
between grid points, so that at the bottom of the ice: 
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(3.18) 



O _ i..(Twfreez ~ To) 

Q cond " Kl (hi/4) 



where k f is the ice thermal conductivity, T Hfreez is the freezing 
point of seawater [=271.3 K] , T 2 is the temperature of the 
second layer of ice, and hj is the ice thickness. 

The storage of latent heat in brine pockets is modeled by 
an internal heat reservoir as in Semtner (1976a) . The 
reservoir accumulates a fraction of the solar radiation that 
penetrates the snow-free ice [Q ( 1-I 0 ) ( l-a jce ) ] . Total heat in 
the reservoir is limited to 50% of the heat needed to melt all 
the ice at that position. This stored heat adds inertia to 
the simulated melting process. The effects of rapid changes 
in thermodynamic forcing are smoothed and the onset of the 
melting season is delayed. A schematic of the fluxes and 
processes incorporated in the thermodynamic portion of the ice 
model is shown in Figure 3.1. The numerical scheme used to 
compute the ice and snow thickness and the temperature at the 
different levels is contained in Semtner (1976a) and will not 
be repeated here. 

The ice dynamics portion of the model uses a simplified 
Kibler (1979) rheology called "bulk viscous" (Washington and 
Parkinson, 1986; Hibler, 1988). This approach was developed 
for use on seasonal time scale ice models with slowly varying 
forcing. Semtner (1987) used this rheology and was able to 
produce annual ice cycles very similar to the Hibler and Bryan 
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Figure 3.1 Schematic of fluxes and ice processes incorporated 
in the thermodynamic portion of the ice model. 
The total heat flux provided by the ocean 
circulation below the mixed layer is shown as Q ocean 
and Q cond is the amount of heat conducted through 
the ice from the surface. The heat reservoir can 
only receive heat when the ice is snow free; 
however, it can provide heat whenever the net heat 
balance at the surface becomes negative. 
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( 1986 ) model which, as noted earlier, used the complete daily 
variable forcing dynamic approach of Hibler ( 1979 ) . The 
smoothed forcing permits an explicit timestepping technique 
to be used vice an iterative approach and the stress tensor 
formulation is much simpler. These two factors reduce the 
computational requirements by approximately 66%. The dynamics 
portion of the ice model predicts changes to the ice thickness 
and concentration from advection, diffusion (a numerical 
requirement) and compression effects. The sea-ice velocity, 
V; , is predicted as follows: 

DVj/Dt = 1/m ( -mfVj + r a + r u - mg Vh + F) (3.18a) 

where D/Dt is the substantial time derivative, f is the 
Coriolis parameter and F is the internal ice stress. r a and 
r w are the wind and ocean current stresses, respectively and 
Vh is the sea surface slope. m is the mass and g is the 
acceleration of gravity. 

The bulk viscous rheology causes the ice to resist 
compression in a plastic fashion but allows unimpeded 
divergence. In the Hibler (1979) notation the stress tensor 
is reduced to: 




(3 . 18b) 



where the strain rate tensor is: 




lftfui , <9ui 



(3 . 18c) 
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The bulk viscosity f is set to zero for divergence (£^>0). 
For convergence, f = P*/2e,j, where P* = h ( x 10 4 [N/m 2 ] . The 
bulk viscosity is not allowed to exceed 2.5 x 10 8 x P* [sec]. 

C. GRID 

The ocean model uses a spherical coordinate system with 
the equator of the system passing through the geographic North 
Pole and the prime meridian situated at 40 degrees East 
longitude. Shifting the poles of the spherical coordinate 
system 90 degrees to the geographic equator, removes the 
singularities of the system and provides more uniform 
horizontal resolution within the Arctic region. The resulting 
horizontal grid interval is approximately 110 kilometers. The 
position of a grid point is defined by latitude, longitude and 
ocean depth. 

The ice portion of the linked model is run on a cartesian 
grid with 110 km spacing. The two grids coincide exactly at 
the pole but are slightly offset at the southern boundaries. 
This introduces a small error in the ice mass balance and salt 
fluxes into the ocean in these areas. Hibler and Bryan (1987) 
determined that the maximum errors introduced into their model 
using a similar scheme were only 10%. A 10% error at the 
southern fringes of the grid was not considered critical to 
the results of this work. 

The gridded domain used in this study was shown in Figure 
1.1. The grid points are staggered on a lattice of type B 
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(Mesinger and Arakawa, 1976), where the temperature, salinity, 
and stream function points are located at the center of a grid 
box. The horizontal components u and v of the velocity are 
at the middle of the vertical edges. The vertical component 
w is computed at the center of the horizontal faces for the 
determination of heat and salt vertical advection, and at the 
box corners for the calculation of u and v (Figure 3.2). 

A total of 13 levels are used in the vertical in order to 
resolve the bathymetric variations shown in Figure 3.3. The 
uppermost level is treated as an isentropic mixed layer of 30 
m. The levels are thinnest near the surface to allow five of 
these levels above the continental shelf, which is at 
approximately 250 m depth. Three more levels are defined down 
to approximately 1100 m depth for representation of the 
Atlantic layer. Table 1 lists the 13 layer thicknesses and 
the depth at the bottom of each layer. 

D. BOUNDARY AND INITIAL CONDITIONS 

The rigid lid approximation requires that w=0 at the upper 
boundary (z=0). The wind stress r[Pa] is specified at the 
surface by 



t( A,<^) — PoA z Qj 2 (u,y) at z = 0 (3.19) 

The small amount of water mass which is transferred at 
the ocean surface is neglected. This includes precipitation 
minus evaporation, ice melt and minor river runoff. However, 
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Figure 3.2 The distribution of levels and grid points. 

Labels to the right and left of the grid points 
refer to ocean and sea-ice variables, 
respectively . 



57 




Figure 3 . 3 



Arctic bathymetry plotted with a 500 m contour 
interval. (from Semtner, 1987). 
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TABLE 1 



VERTICAL LAYER THICKNESSES AND CORRESPONDING DEPTHS 



LAYER 


THICKNESS (M) 


DEPTH (M) 


1 


30 


30 


2 


20 


50 


3 


40 


90 


4 


60 


150 


5 


100 


250 


6 


183 


433 


7 


299 


732 


8 


391 


1123 


9 


487 


1610 


10 


529 


2139 


11 


563 


2702 


12 


585 


3287 


13 


589 


3876 



the resultant changes in salinity from these processes are 
accounted for by applying prescribed negative salt fluxes at 
the surface. The prescribed salt fluxes, as well as the salt 
flux introduced during simulated freezing, are accounted for 
by the salt flux density Q s [K g m ' 2 s* 1 ]: 

Qg = Pwli zff = — /9wSS W (3.20) 

where SSW [m/s] is the sum of all fluxes of fresh water at 
the surface. 

A heat flux density Q 7 [Wm' 2 ] specified at the surface (z=0) 
affects the temperature field in the same way: 

Q T = K z fI=SSH < 3 - 21 ) 
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where SSH [W/m 2 ] is the net energy flux density at the 
surface . 

Vertical motion is forced at the ocean bottom as the 
currents respond to the vertical contours of the bottom slope. 
Bottom friction is determined using a quadratic drag law in 
the same manner as Weatherly (1972). Vertical advective 
fluxes of heat, salt and momentum at the bottom are taken to 
be zero. The vertical diffusive fluxes of heat and salt at 
the bottom are also zero. 

PoK z ^(T,S) = 0 at z = — Ii(A,0) (3.22) 

The horizontal boundaries coincide (to the limits of the 
resolution) with the geographic coastlines and islands of the 
Arctic Ocean and adjacent seas. A no-slip condition, i.e., 
(u,v)=0, is imposed at lateral walls, with no heat or salt 
fluxes across them. 

Each model run simulates the ocean and ice cover for one 
year, producing output fields for the end of each month. The 
decade from 1971-1980 was simulated by conducting ten 
successive runs using the final (end December) ocean and ice 
condition from the previous year’s simulation, as the start 
condition for the current simulation. Since the analysis 
begins with 1971, December 1970 would be used as the very 
first start condition. However, because the model requires 
some time to achieve a steady state after each modification, 
a three year run-up was used. That is, the runs were started 
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at 1968, giving the model a chance to stabilize before the 
data from 1971-1980 were produced and subsequently extracted 
for analysis. 

E. FORCING FIELDS 

This section describes the various forcing fields and data 
used to drive the coupled model. Several formulas are 
presented which contain empirical constants. While the values 
used for these constants represent current estimates, their 
true values remain uncertain. 

Monthly averages of observed atmospheric surface pressure, 
specific humidity, long and short wave radiation and air 
temperature were provided by Dr. J. Walsh. The pressure 
fields were obtained from the NCAR daily analysis of surface 
pressure. This analysis is based on daily reports from 
stations on the periphery of the Arctic Basin and several ice 
stations. Recent comparisons of these analyses with Arctic 
Buoy data show that the average error is limited to 1-2 mb 
(Walsh, private communication) . Temperature fields were also 
obtained from periphery meteorological stations and ice 
stations however they were compiled as monthly averages. The 
temperatures were assumed to be 10 m values. Specific 
humidity at 10 m, q 10 , was calculated using the prescribed 
temperature fields and an assumed relative humidity of 90%. 
Specific humidity at the surface, q s , was calculated using the 
surface temperature produced by the model and assuming 
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saturation. The long and shortwave radiation fields were 
determined in a similar manner to Parkinson and Washington 
( 1979 ) . Radiation was calculated for each hour as a function 
of prescribed zenith angle and cloud cover. The cloud cover 
varied monthly over an annual cycle but was constant across 
the domain. 

Since the prescribed cloud cover was probably the weakest 
of the assumptions made, Walsh has noted (private 
communication) that the long and shortwave radiation fields 
probably contain the most error of all the atmospheric forcing 
data. The specific humidities are limited by the 10 m height 
and 90% relative humidity assumptions, but are still 
considered reasonable. The pressure and temperature fields 
appear quite accurate as confirmed by monitoring stations. 

These were the same forcing data used in Walsh et al., 
(1985). That work, as well as that of Hibler and Walsh 
(1982), showed that inclusion of the dynamic influences from 
interannually varying atmospheric forcing provided improved 
simulations of the interannual variations in the sea-ice 
cover. Time interpolation on a daily basis was conducted 
using a cubic spline method similar to Maykut and Untersteiner 
(1971). The smoothed daily forcing, as opposed to the 
observed daily forcing, was used in order to employ a 
simplified Hibler ice dynamics parameterization as used in 
Semtner (1987). This reduced the computational load of the 
ice portion of the model considerably. 
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Wind stress on the ice and open ocean was computed by 
conventional aerodynamic bulk formula methods in the same 
manner as Semtner (1987) : 

r a = C d/ > a ||V 6 ||V 6 [B] (3.23) 

where C d is the drag coefficient (1.75 x 10' 3 ) . This value is 
larger than the value used by Hibler and Bryan (1987) to 
compensate for the use of time-averaged wind. B is a rotation 
matrix which shifts the forcing vector 25 degrees to the left 
of the surface geostrophic wind when the stress acts on the 
ice. When the wind stress is applied to the open ocean, B was 
set to one (no deflection) . The wind stress on the ice- 
covered ocean was not computed explicitly but was computed 
interactively through the ice (See Section G. Ice-Ocean Model 
Coupling . ) 

Mass inflow through the Faroe-Shetland Channel and the 
Bering Strait and matched outflow through the Denmark Strait 
and the Canadian Archipelago were prescribed. These values 
are the same as in Semtner (1987). The mass and T,S values 
of the inflow, and the mass of the outflow are invariant in 
time, but the T,S properties of the outflow are dependent on 
the simulated fields produced by the model. The inflow of 
fresh water from major rivers around the margins of the Arctic 
Basin is specified on the basis of monthly values tabulated 
by Cattle (1985). Cubic splines are used to interpolate to 
daily values. 
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The flux of sensible heat, Q H (positive downward) , is 
computed by the standard aerodynamic bulk formula method 
(e.g. , Gill, 1982) . 

Qr = Cpj||Vg||(T a — T w ) (3.24) 

where C p is the specific heat of air (1004 J/kg/K) for dry 
air (Huschke, 1959) , C H is a dimensionless transfer 
coefficient for sensible heat called the Stanton Number (1.75 
x 10" 3 (Maykut, 1978), and T w is the water surface temperature. 

The flux of latent heat, Q LE , is computed in a similar 
fashion : 

Q L E = PaLgCgHX sll ~ f lsurf) (3.25) 

where Lg is the latent heat of vaporization or sublimation 
(2.58 x 10 6 J/kg) . C E is a dimensionless transfer coefficient 
called the Dalton Number (1.75 x 10' 3 (Maykut, 1978)), and q 10 
and q surf are the specific humidities at 10 m and at the 
surface, respectively. 

F. SOLUTION 

A comprehensive description of the solution methods for 
the ocean model is contained in Semtner (1986b). The 
following brief summary provides a general overview of the 
procedures . 

The numerical solution of Equations 3.2, 3.3, 3.4, 3.7, 

3.8, 3.9, 3.11, 3.12 and 3.13 is conducted in several steps. 
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The horizontal velocity is decomposed into a baroclinic 
(effect of vertical shear) and a barotropic (vertically 
averaged) component: 



(u,v) = (u\v') + (u,v) (3.26) 

The vertically averaged velocity components (u,v) are written 
in terms of a stream function i> . The vertically integrated 
flow is non-divergent which guarantees the existence of 
The stream function represents the integrated mass transport: 



-If 0 , 1 dv 

U = K J_ h ud2aB -Er;3? 

- i r° , _ i ot 

v = E J-h Vdz = ^ 3,27 

A prediction equation for the stream function, independent of 
the surface pressure is then derived; first by vertically 
averaging and then by taking the curl of the horizontal 
momentum Equations 3 . 2 and 3.3: 
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where G # and represent all the nonlinear and viscous terms 
in Equations 3.2 and 3.3. 

Equation 3.28 essentially relates the time derivative of 
the Laplacian of the stream function to the curl of the 
vertically averaged sum of all the forcing terms in the 
Navier-Stokes equations. The solution requires inversion of 
a second-order differential operator in a closed or multiply 
connected domain. The stream function ip must be a constant 
along coastlines to avoid mass transport across the 
boundaries. \p is held constant in time along the continent 
boundaries. On the islands, \p varies in response to the 
changing circulation in the same manner as Takano (1974), but 
with some modification to allow for variable bottom 
topography. Essentially this involves integrating the curl 
of the vertically averaged momentum equations around each 
island. The condition is then applied so that the line 
integral of the surface pressure gradient, VP S , around each 
island is zero which produces an equation to predict the 
change of ip on each island. 

A prediction equation for the vertical shear of velocity 
is found by differentiating the equation of horizontal motion 
with respect to depth z, and applying the hydrostatic relation 
to eliminate the atmospheric pressure P a . The baroclinic and 
barotropic components of the velocity are then added to obtain 
the total u and v. 
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The prediction equations for temperature and salinity 
(3.11) and (3.12) are much easier to solve. A "leapfrog" 
timestep is used for advection and a "forward" timestep is 
used for diffusion to maintain numerical stability. 
Additionally, a forward timestep is applied every 10 steps 
for the advection calculations. This was an efficient method 
to suppress the computational mode of the leapfrog 
timestepping technique. 

Two different timesteps were used based on the different 
rates of change of temperature, salinity, and velocity in the 
ocean model. Temperature and salinity calculations employed 
6-hour timesteps while the faster varying velocity field used 
a shorter timestep of 22.5 minutes. This reduced the total 
number of calculations required, thereby speeding up the 
numerical integration. The same technique was used 
successfully by Bryan (1984). Different timesteps were also 
used in the ice model. Thermodynamic timesteps were 15 
minutes and dynamic timesteps were 15 seconds to maintain 
numerical stability as in Semtner (1987). 

G. ICE-OCEAN MODEL COUPLING 

The ocean and sea-ice portions of the linked model are 
coupled through momentum, heat and salt exchanges. The 
momentum transfer involved with the coupling occurs at the 
ice-ocean interface. The ice to water momentum transfer is 
modeled by a linear drag formula: 
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(3.29) 



= A 



<9V 



z^r 1 ' C i0 (Vi — V w ) 



where V, is the ice velocity, V u is the ocean velocity in the 
first (mixed) layer and the linear drag coefficient C jo is 
assigned a value of 27.5 x 10' 3 following Semtner (1987). If 
both ice and open water are present in a grid space, the wind 
stress is applied to the open water with Equation 3.19 and the 
ice to ocean stress is computed by Equation 3.29. The average 
stress is then computed by weighting the two contributions 
appropriately by the ice concentration. 

The oceanic heat flux into the mixed layer is provided by 
the ocean model . The computation of the mixed layer 
temperature at each time step and grid point is made in three 
parts : 

- The first computation isolates the first layer 
thermodynamically and dynamically from the rest of the 
ocean. The temperature of the layer T ml , is determined 
simply from the mixed layer/atmospheric heat balance. If 
T ml reaches the freezing point of seawater, sea-ice is 
formed. The change in temperature from the preceding time 
step is equivalent to a heat flux density Q ml (W/m 2 ) at the 
surface : 



Qml = PwCwhml (3.30) 

where h ml is the depth of the mixed layer (30 m) . 

- The second step uses the Q ml calculated in Equation 3.30 
as an upper surface thermodynamic forcing for the 
computation of vertical heat diffusion in the ocean. The 
ocean distributes this input of energy by horizontal 
diffusion, advection and convection. A new temperature 
T dynamic then obtained for the mixed layer. 
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- The final step ensures that energy is conserved. If sea- 
ice is present at the grid point being worked with, and 
if T dynami c is higher than the freezing point of seawater, 
T W freez' then enough ice is melted to decrease T dynamjc back to 
the freezing temperature. In the rare cases where T dyriamjc 
is actually lower than the freezing temperature, ice will 
accrete to the bottom of the existing ice until the latent 
heat of freezing provides enough heat to return the layer 
to the freezing temperature. 

Salt exchanges are computed by prescribing the salt flux 
based on the amount of sea-ice melting and accretion, 
precipitation, evaporation and glacial melt. When no sea-ice 
is present, the mixed layer salinity is computed from Equation 
3.12 and 3.20 by specifying the fresh-water flux at the 
surface, SSW. 

When sea-ice is present, ice accretion and melting produce 
positive and negative salt fluxes, respectively. It is 
assumed that upon freezing, 70% of the initial salt content 
of the seawater is rejected (Semtner, 1987). Therefore the 
salt flux density Q salt (kg s^m' 2 ) at the upper surface (z=0) is 
given by: 



Qsait - Kzff 






melt + 0.70/? -dlAi il i ) 



s ri ow 






-p w [(P-E)(l-Ai)+R] (3.31) 



where A f is the area of the ice cover; = 1 if snow is 
melting and 0 otherwise. P is the precipitation rate, E is 
the evaporation rate and R is the amount of glacial melt. 

Table 2 is a summary of the constants used in the model. 
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TABLE 2 



CONSTANTS 



Time steps ocean T and S 

ocean velocity 
ice T 

ice velocity 



6 hr 

22 . 5 min 
15 min 
15 sec 



Arakawa type B grid 
g , acceleration of gravity 

earth's speed of rotation 
earth’s radius 

vertical eddy viscosity (u,v) 
horizontal eddy viscosity (u,v) 
vertical eddy diffusivity (T,S) 
horizontal eddy diffusivity (T,S,ice) 
density of air 
density of snow 
density of ice 

specific heat capacity of ice and snow 
specific heat capacity of air 
Stanton Number 
Dalton Number 

thermal conductivity of ice 
thermal conductivity of snow 
latent heat of fusion for snow 
latent heat of fusion for ice 
latent heat of evap. (and sublim.) 
melting point for snow 
melting point for snow 
surface emissivity 
Stef an-Bol tzman constant 
f max # maximum limit for bulk viscosity 
C d , air/sea and air/ice drag coefficient 
C jo , ice/water linear drag coefficient 
albedo of ice 
albedo of snow 
albedo of melting snow 
albedo of open ocean 



n 

r e 

a, 

\ 

f>. 

fit 

P\ 

C 

C„ 



W 

h 



9.8 m/s 

7.272 x 10' 5 rad/s 
6.37 x 10® m 
0.3 x 10' 4 m 2 /s 
12 x 10 4 m 2 /s 
0.3 x 10' 4 m 2 /s 
2 x 10 3 m 2 /s 

1 . 3 kg/m 3 
330 kg/m 3 
9 00 kg/m 3 

2.09 kJ/ (kg K) 

1.004 kJ/ (kg K) 

1.75 x lO'^ 3 
1.75 X 10' 3 
4.068 W/(m K) 

0.31 W/fm K) 

110 MJ/m 3 

301 MJ/m 3 
2.58 x 10 6 J/kg 

0 °C 
-0.1 °C 

1 

5 ; 67 x 10'® W/(m 2 K 4 ) 
P* x 2.5 x 10® sec 
1.75 X 10' 3 
27.5 X 10' 3 
0 . 70 
0.80 
0.745 
0 . 07 
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IV. INTERACTIVE OCEAN EXPERIMENT 



This experiment was conducted to determine the importance 
of incorporating a fully interactive prognostic ocean model 
into a linked ice-ocean numerical model. Many numerical 
models of sea ice have been developed which do not include a 
prognostic ocean component. This experiment was designed to 
test the validity of that omission. The importance of 
including a prognostic ocean was judged with respect to the 
accuracy of the simulation of the interannual and annual 
variations of sea-ice cover in the Arctic. Evaluation was 
based on comparisons between observed fields of ice 
concentration and simulated ice concentration fields produced 
by the model . 

A. EXPERIMENTAL SETUP 

The linked ice-ocean model was initially run for ten years 
with a fully interactive ocean and monthly varying atmospheric 
forcing. All oceanic variables were saved on a regular basis 
from this first run, in order to examine their variability and 
to calculate a ten-year average annual cycle. The second run 
used the same atmospheric forcing, integration period and ice 
portion as in the first run. However the interactive, 
prognostic ocean portion was replaced by prescribing the 
average annual cycle from the first run. The outputs from 
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both runs were analyzed to determine each model's skill in 
reproducing the observed annual cycle and interannual 
variations of ice concentration. Modifications to the Semtner 
(1987) model included: 

- Run 1 



Changing the atmospheric forcing from a 30-year 
mean annual cycle to an interannually varying 
cycle of observed monthly averages. 

Defining the initial conditions for each modelled 
year as the end condition for the previous year 
vice averaged or constant fields. This permitted 
consecutive integration of as many years as 
desired. 

Storing and then averaging simulated oceanic heat 
flux and layer 2 ocean velocities at every time 
step for ten years. The oceanic heat flux saved 
for each grid box was the summation of all the 
heat the ocean model provided to the mixed layer 
in that box. Layer 2 velocities (corresponding 
to a depth of 4 0 m) were used because they 
simulated the geostrophic currents quite well. 
The current velocities in the mixed layer were 
unrealistic because all the induced Ekman flow was 
confined within that layer. A new data file 
containing the 10-year average annual cycle of 
these variables was then created. 



- Run 2 



- Incorporating the 10-year average annual cycle of 
oceanic heat flux and layer 2 velocities into the 
model's spline smoothing routine. This permitted 
the continued use of the modified Hibler (1979) 
ice dynamics. 

- Replacing the ocean model with a parameterization 
of the average oceanic forcing (thermodynamic and 
dynamic) . This forcing was applied to the base 
of the mixed layer. 



72 



B. RESULTS 



The average computing time for a single simulation year 
using the fully prognostic, interactive ocean in run 1 was 575 
seconds on the CRAY XMP. The average time for a single year 
in run 2 with the prescribed ocean was 275 seconds. For 
labeling and descriptive purposes, the observed data were 
designated as A data, Run 1 output as B data and Run 2 output 
as C data. 

The simulated ice concentration contours for all 120 
months were produced for B and C in the same manner as for the 
observed fields. As representative examples, Figures 4.1 - 
4.6 show the ice concentration contours for months 76 and 93 
(Apr 77 and Sep 78) for all three data sets (A, B and C) . 
These contours are typical for minimum and maximum ice extent 
periods . 

Layer 2 velocities (Figures 4.7 - 4.10) and contours of 
oceanic heat flux into the mixed layer (Figures 4.11 - 4.14) 
for each month for B and the corresponding average fields for 
C are also shown. Heat units are degrees C per second, 
applied to the volume of the mixed layer within each grid box 
(3.63 x 10 11 m 3 ) . 

Ice thickness was another useful predicted variable. The 
ice thickness distribution provided a good indication of how 
the ice dynamics in combination with the ice rheology modified 
the ice cover within the high ice concentration regions. 
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Figure 4 . 1 



Observed ice concentration (A) in tenths for month 
76 (April, 1977). Contour line labelled -1.5 is 
coastline. 







Figure 4.2 Simulated ice concentration using prognostic ocean 
model (B) for month 76 (April, 1977). Contour 
line labelled -1.6 is coastline. Contours are in 
tenths . 
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Figure 4 . 3 Simulated ice concentration using 10-year mean 
ocean data (C) for month 76 (April, 1977). 
Contour line labelled -1.6 is coastline. Contours 
are in tenths. 
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Figure 4.4 Observed ice concentration (A) in tenths for month 
93 (September, 1979). 
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Figure 4.5 Simulated ice concentration using prognostic ocean 
model (B) for month 93 (September, 1979) . 
Contours are in tenths. 
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Figure 4.6 Simulated ice concentration using 10-year mean 
ocean (C) for month 93 (September, 1979). 
Contours are in tenths. 
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Figure 4.7 



10-year 

currents 



mean, level two (40 m 
for April. Maximum vector 



depth) ocean 
0.246 m/s. 
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Figure 4.8 Simulated level two (40 m depth) ocean currents 
for month 76 (April, 1977) . Maximum vector 0.162 
m/s . 
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Figure 4.9 10-year mean level tv/o (4 0 m depth) ocean currents 
for September. Maximum vector 0.152 m/s. 
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Figure 4.10 Simulated level two (40 m depth) ocean currents 
for month 93 (September, 1979). Maximium vector 
0.344 m/s. 
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Figure 4.11 10-year mean ocean heat flux through the 30 m 
mixed layer for April. Units are degrees 
C/sec/cn 2 scaled by 1 x 10 8 . Conversion factor to 
W/m 2 is 1.25 x 10 s . L and H indicate relative lows 
and highs respectively. 
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Figure 4.12 Simulated ocean heat flux through the 30 m mixed 
layer for month 76. Units are degrees C/sec/crc 2 
scaled by 1 x 10 s . Conversion factor to W/m 2 is 
1.25 x 10 e . L and H indicate relative lows and 
highs .respectively. 
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Figure 4.13 10-year mean ocean heat flux through the 30 m 
mixed layer for Sep. Units are degrees C/sec/cm 
scaled by 1 x 10 e . Conversion factor to W/m 2 is 
1.25 x 10 8 . L and H indicate relative lows and 
highs respectively. 



86 




Figure 4.14 Simulated ocean heat flux through the 30 m mixed 
layer for month 93. Units are degrees C/sec/cm 2 
scaled by 1 x 10 8 . Conversion factor to W/m 2 is 
1.25 x 10 8 . L and H indicate relative lows and 
highs respectively. 
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Representative ice thickness contours for B and C in months 
76 and 93 are shown in Figures 4.15 - 4.18. 

In order to gain a more quantitative appreciation of the 
differences between the various outputs (A, B and C) , a series 
of comparative graphs were produced. The total monthly ice 
area for each of the four regions was determined by 
multiplying the ice concentration at each grid point by the 
grid box area and summing over the region (Figures 4.19 - 
4.22). The observed ice concentration (A) was subtracted from 
the B and C fields to show the amount each simulated field 
differed from the observed (Figures 4.23 - 4.26). The 
absolute sum of the areal differences for B and C in each 
region was then calculated. A measure of the improvement in 
simulating the annual cycle of total ice area in each region 
was then possible. The absolute sum of B-A was consistently 
less than the absolute sum of C-A indicating that the B data 
were closer to the observed data and on average more accurate 
than the C data. The results from these calculations are 
shown in Table 3. 

Figures 4.19 - 4.22 indicated that both B and C had 
obvious biases in total ice area. In order to examine the 
interannual variations, these biases, as well as the mean 
annual cycle of ice cover, were removed. The ten-year average 
annual cycles of A , B and C were calculated (Figures 4.27 - 
4.30) and subtracted from each of their respective data fields 
to produce anomaly (difference from annual cycle) fields of 
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Figure 4.15 Simulated thickness contours using prognostic 
ocean (B) for month 76. Contours are cm of ice 
thickness . 
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Figure 4.16 Simulated thickness contours using 10-year mean 
ocean (C) for month 76. Contours are cm of ice 
thickness . 
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Figure 4.17 Simulated thickness contours using prognostic 
ocean (B) for month 93. Contours are cm of ice 
thickness . 
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Figure 4 . 18 Simulated thickness contours using 10-year mean 
ocean (C) for month 93. Contours are cm of ice 
thickness . 
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A1/B1/C1 




Figure 4.19 Time series of total ice area in region 1. 

Observed (A) is the solid line, prognostic 
ocean model (B) is the dashed line and 10- 
year mean ocean model (C) is the dash-dot 
line. X axis is months from Jan 1971 to Dec 
1980. Y axis is 10's of km 2 . 
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A2/B2/C2 




Figure 4.20 Tine series of total ice area in region 2. 

Observed (A) is the solid line, prognostic 
ocean model (B) is the dashed line and 10- 
year mean ocean model (C) is the dash-dot 
line. X axis is months from Jan 1971 to Dec 
1980. Y axis is 10's of km 2 . 
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A3/B3/C3 




Figure 4.21 Time series of total ice area in region 3. 

Observed (A) is the solid line, prognostic 
ocean model (B) is the dashed line and 10- 
year mean ocean model (C) is the dash-dot 
line. X axis is months from Jan 1971 to Dec 
1980. Y axis is 10's of km 2 . 
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A4/B4/C4 




Figure 4.22 Time series of total ice area in region 4. 

Observed (A) is the solid line, prognostic 
ocean model (B) is the dashed line and 10- 
year mean ocean model (C) is the dash-dot 
line. X axis is months from Jan 1971 to Dec 
1980. Y axis is 10's of km 2 . 
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BDIF1/CDIF1 




Figure 4.23 Time series of simulated minus observed total area 
in region 1. B-A is the solid line, C-A is the 

dashed line. X axis is months from Jan 1971 to 

Dec 1980. Y axis is 10's of km 2 . 
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BDIF2/CDIF2 




Figure 4.24 Time series of simulated minus observed total area 
in region 2. B-A is the solid line, C-A is the 
dashed line. X axis is months from Jan 1971 to 
Dec 1580. Y axis is 10 ’s of km 2 . 
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BDIF3/CDIF3 




Figure 4.25 Time series of simulated minus observed total area 
in region 3. B-A is the solid line, C-A is the 

dashed line. X axis is months from Jan 1971 to 

Dec 1980. Y axis is 10's of km 2 . 



99 



BDIF4/CDIF4 




Figure 4.26 Time series of simulated minus observed total area 
in region 4. B-A is the solid line, C-A is the 
dashed line. X axis is months from Jan 1971 to 
Dec 1980. Y axis is 10's of km 2 . 
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A1MN/B1MN/C1MN 




Figure 4.27 10-year mean annual cycle of total ice area in 
region 1. Observed (A) is the solid line, 

simulated (B) is the dashed line and simulated (C) 
is the dash-dot line. X axis is months and Y axis 
is lo's of km 2 . 
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A2MN/B2MN/C2MN 




Figure 4.28 10-year mean annual cycle of total ice area in 
region 2. Observed (A) is the solid line, 

simulated (B) is the dashed line and simulated (C) 
is the dash-dot line. X axis is months and Y axis 
is 10 ' s of km 2 . 
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A3MN/B3MN/C3MN 




Figure 4.29 10-year mean annual cycle of total ice area in 
region 3. Observed (A) is the solid line, 
simulated (B) is the dashed line and simulated (C) 
is the dash-dot line. X axis is months and Y axis 
is 10' s of km. 
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A4MN/B4MN/C4MN 




Figure 4.30 10-year mean annual cycle of total ice area in 
region 4. Observed (A) is the solid line, 

simulated (B) is the dashed line and simulated (C) 
is the dash-dot line. X axis is months and Y axis 
is 10' s of km 2 . 
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TABLE 3 



DIFFERENCE BETWEEN THE ABSOLUTE SUM OF SIMULATED 
AND OBSERVED ICE AREA 



REGION 


B - A 


0 

1 

> 


CHANGE 


% IMPROVEMENT 


1 


.4391 


.4472 


.0081 


2% 


2 


.1645 


. 1729 


.0084 


5% 


3 


. 2945 


.3131 


.0186 


6% 


4 


.2043 


.2935 


.0892 


43% 



area in km 2 x 1 x 10 8 

ice concentration for A,B and C in each region. The B and C 
anomaly time series are shown with the A anomaly time series 
for each region in Figures 4.31 - 4.38. The ability of the 
B and C anomaly series to follow the A series was a measure 
of their ability to model interannual variability. A 
quantitative appreciation of this ability was made in a 
similar manner as with the annual cycle. The A anomaly field 
was subtracted from both B and C and the absolute sum of each 
of the anomaly difference fields was then calculated. The 
results are shown in Table 4. Correlations were also 
calculated between the observed and modelled ice area fields 
for both the annual cycle and anomaly time series. The 
calculated correlation coefficients are shown in Table 5. 

Several hundred contour plots and graphs were examined and 
compared including observed and simulated ice concentration 
and ice thickness, simulated ocean heat flux, simulated layer 
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A1DIFM/B1DIFM 




Figure 4.31 Difference from mean annual cycle or anomaly time 
series of total ice area for region 1. Observed 
anomaly is the solid line, simulated (B) anomaly 
is the dashed line. X axis is months, Y axis is 
10's of km 2 . 
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A1DIFM/C1DIFM 




Figure 4.32 Difference from mean annual cycle or anomaly time 
series of total ice area for region 1. Observed 
anomaly is the solid line, simulated (C) anomaly 
is the dashed line. X axis is months, Y axis is 
10's of km 2 . 
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A2DIFM/B2DIFM 




Figure 4.33 Difference from mean annual cycle or anomaly time 
series of total ice area for region 2. Observed 
anomaly is the solid line, simulated (B) anomaly 
is the dashed line. X axis is months, Y axis is 
1 0 1 s of km 2 . 
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A2DIFM/C2DIFM 




Figure 4.34 Difference from mean annual cycle or anomaly time 
series of total ice area for region 2. Observed 
anomaly is the solid line, simulated (C) anomaly 
is the dashed line. X axis is months, Y axis is 
1 0 ' s of km 2 . 
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A3DIFM/B3DIFM 




Figure 4.35 Difference from mean annual cycle or anomaly time 
series of total ice area for region 3. Observed 
anomaly is the solid line, simulated (B) anomaly 
is the dashed line. X axis is months, Y axis is 
10 ' s of km 2 . 
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A3DIFM/C3DIFM 




Figure 4.36 Difference from mean annual cycle or anomaly time 
series of total ice area for region 3 . Observed 
anomaly is the solid line, simulated (C) anomaly 
is the dashed line. X axis is months, Y axis is 
10's of km 2 . 
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A4DIFM/B4DIFM 




Figure 4.37 Difference from mean annual cycle or anomaly time 
series of total ice area for region 4 . Observed 
anomaly is the solid line, simulated (B) anomaly 
is the dashed line. X axis is months, Y axis is 
1 0 1 s of km 2 . 
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A4DIFM/C4DIFM 




Figure 4.38 Difference from mean annual cycle or anomaly time 
series of total ice area for region 4 . Observed 
anomaly is the solid line, simulated (C) anomaly 
is the dashed line. X axis is months, Y axis is 
10's of km 2 . 
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TABLE 4 



DIFFERENCE BETWEEN ABSOLUTE SUM OF SIMULATED AND OBSERVED 

ICE AREA ANOMALY 
(MEAN ANNUAL CYCLE REMOVED) 



REGION 


B - A 


C - A 


CHANGE 


% IMPROVEMENT 


1 


.1476 


. 1818 


. 0342 


23% 


2 


. 0864 


. 0985 


. 0121 


14 % 


3 


. 1355 


. 1689 


. 0334 


25% 


4 


. 0888 


. 1157 


. 0268 


30% 



area in km 2 x 1 x 10 E 



TABLE 5 

CORRELATIONS OF ICE AREA BETWEEN SIMULATION B AND OBSERVED 
(CORR B) AND SIMULATION C AND OBSERVED (CORR C) 



ICE AREA TIME SERIES 



REG 


STDEV OBS 


STDEV B 


STDEV C 


CORR B 


CORR C 


1 


. 3169 


. 7776 


.8246 


.9126 


.8852 


2 


. 1434 


. 3729 


. 3996 


.9271 


. 8807 


3 


.4283 


. 6755 


. 6910 


.9461 


.9106 


4 


.2292 


. 3667 


.4220 


.9186 


.8332 



ANOMALY TIME SERIES 



REG 


STDEV OBS 


STDEV B 


STDEV C 


CORR B 


CORR C 


1 


. 1146 


. 2220 


.2701 


. 4936 


.3963 


2 


. 0590 


. 1660 


. 1821 


.7761 


.5981 


3 


. 1376 


. 1518 


. 1677 


. 5375 


. 3039 


4 


.0994 


. 1311 


. 1418 


.6172 


. 3332 



Standard deviation values x 1 x 10 6 

B simulation was with fully interactive ocean model 
C simulation was with 10-year mean cycle ocean 
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2 velocities and the time series' described above. On 
conclusion of this review, the following general observations 
were made: 

- The region of transition from 10/10 to 0/10 ice 
concentration was narrower in both B and C than in A, 
i.e. , the ice edge boundary was sharper, more well defined 
in B and C than in A. 

- Both B and C overpredicted the ice coverage in winter, 
especially in the eastern Barents Sea (contained in region 
4) . B was usually closer to the observed coverage than 
C. 



- Both B and C underpredicted the coverage in summer, 
especially in the Kara Sea to the east and the Beaufort, 
Chukchi and East Siberian Seas to the west. B was usually 
closer to the observed coverage than C in the eastern 
Arctic but the improvement in the central and western 
Arctic was minimal. 

- The general shape of the ice edge was reasonably well 
simulated by both B and C. The Kara and Barents Seas were 
notable exceptions. 

- The simulated surface (level 2) currents matched the 
observed annual mean currents (Figure 2.4) quite well. 
All the major currents were represented and their 
velocities appeared reasonable; however, some errors were 
evident in the secondary currents. The West Spitsbergen 
Current was not represented and the Irminger Current was 
flowing opposite to the observed mean direction. There 
were significant differences between the monthly varying 
fields and the 10-year average cycle fields of level 2 
ocean velocities. The fields contained the same current 
features however these features differed in both strength 
and position. 

- The ocean heat flux into the mixed layer was highly 
variable between regions and changed dramatically with 
the seasons. Winter values were highest in the Norwegian 
and Greenland Seas where it averaged approximately 180 
W/m 2 . Exceptional cases of heat flux in excess of 375 W/m 2 
were observed occasionally in a localized area just south 
of Svalbard. Spring and summer heat fluxes in the 
Norwegian and Greenland Seas were generally weakly 
positive or negative. Isolated summer heat flux values 
less than -180 W/m 2 were observed just north of Iceland in 
some years. A similar positive to negative ocean heat 
flux cycle was observed in the western arctic although it 
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was generally confined to within 400 km of the Alaskan 
shoreline. The simulated ocean heat flux under the 
central pack in the Canadian and Eurasian Basins was 
negligible throughout the year. Similar to the 
velocities, there were significant differences between 
the monthly varying fields and the 10-year average cycle 
fields in both the strength and position of the main 
features . 



The simulated pack ice thicknesses in the central Arctic 
Basin were 1-2 m or 50-100% less than the average observed 
values, particularly in summer. B thicknesses tended to 
be larger than those for C. The simulated thickness data 
indicated that some ridging was occurring in the models 
poleward of Greenland and Ellesmere Island as expected, 
however the thickness values in this region were too low 
by a factor of three or more. The observed thickness 
fields used for comparison are shown in Figures 4.39 - 
4.42. 

Region 1 and 2 had maximum ice areas limited by their 
boundaries. Region 4's minimum ice area was limited by 
zero (it becomes totally ice free). Region 3 came close 
to being ice free especially from 1978-1980. These 
minimum and maximum limits reduced the variability of the 
total ice area; however considerable yearly variability 
was still evident. Annual differences in ice areas in the 
observed data varied from approximately 15% of the mean 
in region 1 and 2, to 40% in region 3, to 100% in region 
4. The simulated data had greater variance with 
interannual ice area differences ranging from 
approximately 40% of the mean in regions 1 and 3, to 60% 
in region 2, to 100% in region 4. 

The mean annual cycles of ice area for A,B and C showed 
that the melt/freeze cycle in the modelled cases was 
approximately in phase with the observed annual cycle but 
of higher amplitude. C generally displayed the largest 
amplitude . 

The time series' of ice area anomalies (difference from 
annual cycle) for both observed and simulated data varied 
significantly from region to region. Some similarities 
could be found in adjacent regions. However, if they were 
not adjacent, the correlations appeared to be very low. 

The degree of monthly variability or "noisiness" in the 
anomaly series increased markedly in regions 3 and 4 
compared to regions 1 and 2. 
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Figure 4.39 Mean ice draft (m) in winter derived from 
submarine upward looking echo sounder data. Draft 
is averaged over three months, centered on the 
first of March (from Bourke and Garret, 1987). 
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Figure 4.40 Mean ice draft (m) in spring derived from 
submarine upv.-ard looking echo sounder data. Draft 
is averaged over three months, centered on the 
first cf June (from Bourke and Garret, 1987). 
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Figure 4.41 Mean ice draft (m) in summer derived from 
submarine upward looking echo sounder data. Draft 
is averaged over three months, centered on the 
first of September (from Bourke and Garret, 1987) . 
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Figure 4.42 Mean ice draft (m) in fall derived from submarine 
upward looking echo sounder data. Draft is 
averaged over three months, centered on the first 
of December (from Bourke and Garret, 1987). 
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Figure 4.43 10-year mean ocean heat flux through the 30 m 
mixed layer for July. Units are degrees C/sec^cm 2 
scaled by 1 x 10 9 . Conversion factor to W/m is 
.125 x 10 9 . L and H indicate relative lows and 
highs respectively. 



121 



The B anomaly time series was a better match to the 
observed (A) anomaly time series than C. There were time 
periods when both B and C were significantly different 
from A but with the same sign. There are also times when 
they both differed from the A anomalies but with opposite 
signs. No clear pattern was evident of a bias to either 
sign . 

Considering the entire 120 month time period, model runs 
using the interactive ocean (B) matched the observed data 
better than the average ocean case (C) in all regions. 
This was true for both total ice area and interannual 
variability. Both the absolute sums of total differential 
area and the correlation coefficients for each region were 
better for B than for C, particularly in region 4. 
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C. DISCUSSION 



The excessive melting which occured in regions 1 and 2 was 
a common feature of both the prescribed ocean model and the 
prognostic ocean model. This suggests that an error or 
weakness common to both model formulations may be responsible. 
One possible weakness is the quite simple representation of 
the mixing processes within the surface layer and the rather 
large vertical resolution of the upper 500 m. Contours of 
ocean heat flux (e.g., Figure 4.43) show large positive values 
in summer when the melting rate is highest. The large amount 
of melting in this season would be expected to produce a 
surface cap of fresh water. The resultant stratification 
would insulate the mixed layer from the heat below and reduce 
this flux. It appears that this mechanism is not being 
represented properly as excessive heat continues to enter the 
mixed layer until late summer. Alternatively, the large 
lateral diffusion coefficients which are necessary for 
numerical stability may be mixing far too much heat throughout 
the intermediate layer. In this case, even if the vertical 
advection is reduced by stratification, there may still be 
enough heat available to cause the excessive summer ice 
retreat . 

A second possibility becomes evident after sequentially 
viewing the ocean heat flux contours over several months. A 
close examination of those areas which appear most in error 
indicates several "hot" spots in winter and spring which 
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correspond to some of the locations of prescribed inflow. 
These same "hot" spots become "cold" spots during summer and 
fall. The constant prescribed inflow temperatures at the 
Bering Strait and Faeroe-Shetland Channel contrast with the 
simulated ocean temperatures produced by the model. The 
simulated temperatures vary monthly following the annual cycle 
of SST resulting in a cycle of anomalous temperatures at the 
inflow positions. Although the temperature of the river 
inflow is set to match the ocean temperature, similar "hot" 
spots occur at the river mouths too, particularly for the 
Mackenzie, Kolyma and Lena Rivers. The most probable causes 
for this appear to be geostrophic adjustment to the negative 
salt fluxes and/or bottom topography interaction. The heat 
from these hot spots is advected by the currents, spreading 
their effect and increasing the spring melting over a 
considerable area. The western boundary areas appear to be 
a good example of this. 

The simulated winter ice cover, although not an exact 
match to the observed, is still a closer match than the summer 
condition. This would seem reasonable because the large 
freezing rates and surface cooling in the winter induce deep 
vertical mixing. This process can be better represented by 
the simple density mixing scheme used in this model than the 
shallow stratification process which occurs in summer. 
Realistic amounts of heat are brought to the surface and the 
simulated ice edge matches the observed quite well. The 
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improvement provided by the prognostic ocean in regions 3 and 
4 is consistent because, if we consider that the winter 
vertical convection is reasonably well accounted for, an 
improvement in the interannually varying currents and ocean 
heat flux should have a direct improvement on the ice cover 
simulation . 

D. CONCLUSIONS 

Inclusion of the interactive ocean produced an improved 
accuracy of simulation of the annual cycle of total ice area, 
particularly in region 4. The model which included the 
interactive, prognostic ocean represented the variable ocean 
forcing (heat and momentum) much better than the average ocean 
model. This suggests that variable ocean forcing is important 
to the evolution of the ice cover in the eastern Arctic. By 
contrast, in regions 1 and 2, the interactive model did not 
improve the representation of the total ice area 
significantly. In these regions, it may be the atmospheric 
forcing, the ice rheology or as discussed above, mechanisms 
controlling mixing of heat in the upper ocean layers, that 
dominate the evolution of total ice area. These are 
parameters which were not changed between experiments. 
Alternatively, there could be errors in the ocean forcing 
which produces a bias in the ice field in both the variable 
ocean and the average ocean. 
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Inclusion of the interactive model produced significant 
improvements in simulating the interannually varying ice field 
in all regions of the Arctic. This indicates that the 
interannual changes in the ocean forcing are a major influence 
on the interannual variability of the ice field over the 
entire Arctic basin. The importance of including an 
interactive ocean in seasonal ice prediction models is 
therefore demonstrated. 

Simulated ice thickness has not been improved 
substantially based on the data presented by Bourke and Garret 
(1987). Again, as for total ice area in the western Arctic, 
modifications to the way atmospheric forcing is handled (e.g., 
albedo, surface drag) or to the ice rheology may be required 
to improve this feature of the ice cover. 
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V. 



ICE STRENGTH EXPERIMENT 



One of the most obvious inaccuracies in the simulated ice 
cover produced by the linked model was the ice thickness as 
discussed in Chapter IV. In the polar pack region, simulated 
thicknesses were approximately one half the observed values. 
Thicknesses in the ice convergence regions north of Greenland 
and the Canadian archipelago were too small by a factor of 
three or more. 

Heavy ice buildup in the convergence regions is the result 
of the large amount of ridging and rafting which occurs there. 
The poor representation of this feature indicated that the ice 
rheology used was not portraying these mechanisms properly. 
Hibler (1979, 1980) was able to simulate reasonable ice 
thicknesses in the strong convergence areas along the north 
coast of Ellesmere Island, but noted that the simulated ice 
buildup in this region and overall spatial ice thickness 
variations were dependent on the ice strength value used. In 
the follow-on work by Hibler and Bryan (1987) , a greater ice 
strength was used to provide more accurate ice velocities at 
the larger grid scale. However, the ice thickness in the 
convergence regions was considerably reduced, appearing quite 
similar to the initial results found in this work. 

Examination of the ice strength value used in the model 
here indicated that the ice strength parameter, P*, was 
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relatively large. This made the ice pack quite stiff, 
reducing the amount of ridging in areas of ice convergence. 
The ice rheology used in this model was the same as in the 
Semtner (1987) model. However, it should be noted that the 
strength parameter actually used by Semtner (1987) differed 
from the description contained in the Appendix to that 
reference. The P* actually used was: 

P* = h x 3 x 10 5 (N/m 2 ) 

where h is the ice thickness. This was 30 times larger than 
the value noted in the appendix. Semtner indicated (private 
communication) that this larger value was required to prevent 
the numerical instability which had occurred in earlier runs 
using a smaller P* value. 

This experiment investigated the effects of reducing P*. 
The objective was to determine the best value for producing 
an accurate simulation of the ice thickness and thickness 
distribution . 

A. EXPERIMENTAL SETUP 

The fully prognostic ice-ocean model (B) was used for this 
experiment. Initial runs with varying P* values indicated 
that changes to the simulated ice fields were dramatic. 
Indeed, changes were evident after only a single month of 
integration. Therefore, single year integrations were 
considered sufficient for the sensitivity study vice a ten 
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year integration period. The year 1974 was selected for the 
runs as it was a fairly average ice cover year. 

Seven runs were conducted. P* values were set as follows: 

- P* = h x 3 x 10 5 

- P* = h x 1 x 10 5 

- P* = h x 5 x 10 4 

- P* = h x 3 x 10 4 

- P* = h x 1 x 10 4 

- P* = h x 5 x 10 3 

- P* = h x 1 x 10 3 

Once an optimum P* was determined, the full ten year 
integration was conducted to produce ice concentration fields 
for comparison with the observed data. 

B. RESULTS 

Numerical instability occurred in this model at the 
reduced P* values similar to Semtner's work (private 
communication) . The instability originated at a grid point 
in the EGC, adjacent to where the grid representation of the 
Greenland coastline was very jagged. The jagged coastline was 
a result of the model's limited resolution. Ocean current 
velocities grew exponentially, producing artificially high ice 
velocities and excessive kinetic energy. The temperature and 
salinity values also became unreasonable. Reduction of the 
model's numerical time steps by a factor of two did not solve 
the problem. Numerical stability was eventually obtained by 
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imposing velocity, temperature and salinity limits. The 
limitations on these prognostic variables were as follows: 
-75.0 < UICE , VICE < +75.0 (cm/sec) 

-2.0 < TSURF < 15.0 (deg centigrade) 

0 < SALINITY < 50 (ppt) . 

Figures 5.1 - 5.6 show the simulated ice thickness 

contours for August 1974, using the P* values from runs 2 
through 7, respectively. In general, as P* was reduced, the 
entire ice pack tended to shift toward the Fram Strait and 
thicknesses within the pack increased. The convergence along 
the north coast of Greenland became more intense eventually 
producing a reasonable simulation of the ridging and rafting 
processes. However the simulated region of maximum 

thicknesses was still thinner and was located further poleward 
than the average thickness fields in Bourke and Garrett (1987) 
(See Figures 4.39 - 4.42). The tendency for the model to 
produce too much meltback of the ice edge in summer was 
exacerbated, even with the increased ice thicknesses. This 
was especially evident in the Bering and East Siberian Seas. 

The thickness gradients became increasingly stronger, 
irregular and spotty as the ice strength was set to the lower 
values. The two lowest P* values produced thickness contours 
which appeared too "noisy" to be realistic. 

A qualitative assessment of all the thickness contour 
plots and comparison with the Bourke and Garrett (1987) data 
indicated that P* = h x 1 x lO 4 provided the most reasonable 
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Figure 5.1 Simulated ice thickness contours (cm) for August, 
1974 . P* = h x 1 x 10 5 . 
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Figure 5.2 



Simulated ice thickness contours (cm) for August, 
1974 . P* = h x 5 x 10 4 . 
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Figure 5.3 Simulated ice thickness contours (cm) for August, 
1974 . P* = h x 3 x 10 4 . 
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Figure 5.4 Simulated ice thickness contours (cm) for August, 
1974. P* = h x 1 x 10 4 . 
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Figure 5.5 Simulated ice thickness contours (cm) for August, 
1974 . P* = h x 5 x 10 3 . 
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Figure 5.6 Simulated ice thickness contours (cm) for August, 
1974 . P* = h x 1 x 10 3 . 
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simulated overall ice thickness and thickness distributions. 
This value is somewhat larger than that used in Hibler 
(1979,1980) (P* = f(h) x 5 x 10 3 ) and smaller than that used 
in Hibler and Bryan (1987) (P* = f(h) x 2.75 x 10 4 ) . Both of 
those were chosen by matching simulated drift rates with 
observed data. The reason for a different optimum value for 
P* is probably due to the higher spatial resolution and 
monthly averaged daily atmospheric forcing used in this model 
vice the observed daily forcing used in the Hibler models. 

The ten-year integration was run with the optimum P* and 
a fully interactive ocean to produce ten years of simulated 
ice concentration fields. Correlations with the observed data 
for both the annual cycle and interannual variations for all 
four regions were calculated. The correlations using the 
reduced P* model were considerably lower than when the 
experiment 1 model with a large P* value was used (Table 6) . 

C. DISCUSSION 

The linked ice-ocean model with a prognostic ocean 
displayed notable sensitivity to changes in the ice strength 
factor of the ice rheology. Reducing the P* value by a factor 
of 30 improved the simulated ice thickness and distribution 
significantly. However, further reductions in ice strength 
appeared to introduce excessive "noise" in these fields. The 
contour lines became very rough and erratic with numerous 
"bullseyes" or spot irregularities. This was probably due to 
the diminished ability of the ice to sustain internal stress. 
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TABLE 6 



CORRELATIONS OF ICE AREA BETWEEN SIMULATION B AND OBSERVED 
(CORR B) AND SIMULATION D AND OBSERVED (CORR D) 



AREA TIME SERIES 



REG 


STDEV D 


CORR B 


CORR C 


CORR D 


1 


.8160 


.9126 


. 8852 


.8610 


2 


.3882 


.9271 


. 8807 


.8866 


3 


. 6709 


.9461 


.9106 


.9322 


4 


.4394 


.9186 


.8332 


.8623 



ANOMALY TIME SERIIES 



REG 


STDEV D 


CORR B 


CORR C 


CORR D 


1 


. 1697 


.4936 


. 3963 


.4294 


2 


. 1499 


. 7761 


.5981 


. 6292 


3 


. 1748 


. 5375 


.3039 


.4536 


4 


. 1788 


. 6172 


. 3332 


.4354 



Standard deviation values x 1 x 10 6 



B simulation was with 
C simulation was with 
D simulation was with 



fully interactive ocean model 
10-year mean cycle ocean 
the reduced P* 
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As the ice strength was reduced, the ice approached a free 
drift condition with little resistance to either shearing or 
compressive forces. In this condition, the ice could respond 
rapidly and easily to any short term or small scale (a few 
hundred km) dynamic forcing resulting in the irregular 
appearance of the thickness contours at the lowest P* values. 
Since the atmospheric forcing was smoothed, these shorter or 
small scale forcing events were probably from the ocean. 

An explanation for the increased ice thicknesses when a 
moderate reduction in ice strength was introduced is that by 
reducing the ice strength, the long term dynamic influences 
of the Transpolar Drift Stream and atmosphere were able to 
push the ice further and faster to the southeast. The 
resultant ice divergence from the Bering and adjacent seas 
created large expanses of open water. Large heat losses to 
the atmosphere associated with this open water and cold air 
temperatures in fall and winter encouraged rapid ice 
production. The large volumes of new ice were then advected 
to the southeast where at the lower strength levels they were 
able to compress and undergo the appropriate ridging and 
rafting processes. This produced the large areas of thin ice 
on the Siberian shelves, increased the average ice pack 
thicknesses and resulted in the stronger thickness gradients 
evident in the simulated ice thickness fields. 

The ice thickness and distribution fields from the reduced 

* 

P model appeared to be so much improved, that the 
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correlations of ice concentration with the observed fields 
were expected to improve also. The subsequent determination 
that these correlations were actually degraded was surprising. 
Ice thickness and distribution appeared to improve 
substantially within the ice pack. However the amount of 
ocean area which experienced a total loss of ice actually 
increased, particularly in the Bering and East Siberian Seas. 
Since this exaggerated the already excessive ice edge retreat 
in summer, and no compensatory improvement in the ice edge 
position was produced in the other seasons, the correlations 
with the observed data decreased. 

The increased ice edge retreat and areas of open water can 
be explained as follows. The large heat loss over the open 
water and the strong positive salt flux from the large 
freezing rate would promote increased convective overturning 
of the ocean in the Siberian shelf region. Strong convection 
would break through the halocline which serves to insulate the 
surface waters from the warmer waters below. Oceanic heat 
flux to the ocean surface would increase, slowing the freezing 
rate as winter progressed and keeping the ice cover relatively 
thin in those regions. Once the increased solar flux in 
spring produced a net positive heat flux at the ocean surface, 
the thin ice would rapidly melt and extensive areas of open 
water would form. 

Hibler and Bryan (1987) describe the complex 
interactions and feedback mechanisms possible, particularly 
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in the marginal ice zone, when the effects of dynamics and an 
interactive ocean are included. They note that a fine balance 
of ice processes appears to exist. The initial freezing 
causes intense overturning. The resultant deeply mixed 
surface waters can withstand considerable amounts of melting 
at a later time without becoming stratified. This description 
is consistent with the process suggested here. 

D. CONCLUSION 

An appropriate function to use for ice strength in the ice 
rheology of this model appears to be: 

P* = h x 1 x 10 4 

The ice thickness and distribution fields produced using this 
function were considerably more realistic in comparison with 
available observations than those produced using a greater ice 
strength. The simulated average pack thickness was doubled 
and strong ice buildup occurred in a manner and in areas 
consistent with observed data. A possible exception was the 
region immediately poleward of Greenland and the Canadian 
Archipelago. Ice thicknesses there were reduced, apparently 
due to advection because there did not appear to be 
significant ocean heat flux. Smaller ice strength values 
produced thickness fields which were considered too noisy. 

Both the annual cycle and interannual variations of ice 
concentration produced using the reduced ice strength were 
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degraded, primarily because the total ocean area which became 
ice free in summer increased from the initial model runs. It 
is proposed that the reason for the increased area of open 
ocean in summer is an increase in the ocean heat flux. This 
occurs as a result of the strong convective penetration of the 
halocline resulting from greater heat losses to the atmosphere 
from regions where the ice has been thinned or removed and 
increased salt fluxes into the ocean caused by the larger 
freezing rates in the same regions. 

These results indicate that the simulated heat flux 
provided by the ocean appears to be a dominant factor 
controlling the differences between simulated and observed ice 
edge position. Inclusion of the interactive ocean in the 
linked ice-ocean model provided a much more realistic 
representation of this flux. Consequently, simulated ice 
fields were much improved over the model which used only an 
average ocean condition. However, the inaccuracies that 
remain in the representation of ocean heat flux appear to be 
very important. 
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VI 



ALBEDO EXPERIMENT 



The excessive meltback of the ice edge in summer was a 
common feature of previous model runs in this work. Similar 
results have been reported by other ice modellers as well 
(e.g., Semtner, 1987; Hibler and Walsh, 1982). The ice 
strength experiment indicated that this may have been due to 
excessive simulated oceanic heat flux in the areas of severe 
melting. Alternatively, excessive atmospheric heat flux into 
the ice cover could also induce an exaggerated meltback of the 
ice edge. One of the most important elements determining the 
atmospheric heat flux into the ice is surface albedo. 

Shine and Henderson-Sellers (1985) noted the differences 
in albedo representations used by several authors (e.g., 
Parkinson and Washington, 1979; Hibler and Walsh, 1982; Manabe 
and Stouffer, 1980) . They showed that a thermodynamic ice 
model very similar to the one incorporated into this model was 
quite sensitive to changes in the surface albedo. Large 
differences between various model thickness predictions could 
be explained by these albedo differences. 

The Shine and Henderson-Sellers (1985) study used a model 
which did not include ice dynamics nor an interactive ocean. 
As noted in Chapter V, inclusion of these in a linked model 
permits a much more realistic representation of the complex 
mixing and feedback processes which occur during freezing and 
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melting. The response of the more elaborate model used in 
this study to surface albedo changes was unknown and therefore 
was considered worthwhile for further investigation. 

A. EXPERIMENTAL SETUP 

The represention of surface albedo used in this model was 
described in Chapter III. Some allowance was made for the 
differences in albedo between ice and snow and for the changes 
which occur when melting commences (see Chapter III) ; however, 
the representation was very simple. 

Ross and Walsh (1987) used a more elaborate surface albedo 
representation in conjunction with a dynamic-thermodynamic 
model to simulate surface albedo in the Arctic. The albedo 
representation was a linear function of air temperature which 
also accounted for the differences between ice and snow. The 
simulated albedos produced with that model matched the 
satellite-derived estimates by Robinson et al., (1986) very 
well. In view of the importance placed on an accurate 
representation of the surface albedo by Shine and Henderson- 
Sellers (1985), the albedo representation of Ross and Walsh 
(1987) was incorporated into the model used in this work as 



follows : 






a snow — 


0.80 


T sfc < “ 5 -° de< 3 C 


a snow — 


0.65 - 0.03 (T sfc ) 


-5.0 < T sfc < 0.0 d 


a snow — 


0.65 


T S fc = 0,0 de 9 c 


- a ice = 


0.65 


T sfc < 0.0 deg C 



144 



a: 



ice 



ice 



0.65 - 0.04 (T ajr ) 



0.45 



0.0 < T ajr < 5.0 deg C 
T air > 5.0 deg C 



a, 



water 



0.10 



where "a" is the albedo value, T sfc is the surface temperature 
of the ice or snow and T air is the air temperature above the 
ice or snow surface. The fraction of solar radiation absorbed 
internally by the ice was 0.35 as recommended by Ross and 
Walsh (1987). 

A ten year integration was conducted using the new albedo 
model and radiation absorption constant. The resultant ice 
concentration fields were again correlated with the observed 
data as in the previous experiments. The ten-year 
correlations of simulated annual ice cycle and interannual 
variations are shown in Table 7. The correlation coefficients 
improved nominally in regions 2 and 4 and worsened nominally 
in regions 1 and 3. 

Comparison of simulated ice concentration contours and ice 
thickness contours from model runs before and after the albedo 
changes showed only trivial differences. It appeared that the 
full model with ice dynamics and an interactive ocean was 
relatively insensitive to changes in albedo. This finding was 
a clear contrast to previous work using simpler ice models 
(e.g., Shine and Henderson-Sellers , 1985; Hibler, 1980). 



Vote that this is a correction to Ross and Walsh (1987) 
to make the albedo a linear function of temperature as 
described in the text. 
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TABLE 7 



CORRELATIONS OF ICE AREA BETWEEN SIMULATION B AND OBSERVED 
(CORR B) AND SIMULATION D AND OBSERVED (CORR D) 



ICE AREA TIME SERIES 



REG 


STDEV E 


CORR B 


CORR D 


CORR E 


1 


.8066 


.9126 


.8610 


.8606 


2 


. 3778 


.9271 


.8866 


.8870 


3 


. 6647 


.9461 


.9322 


.9321 


4 


.4387 


.9186 


.8623 


.8648 



ANOMALY TIME SERIES 



REG 


STDEV E 


CORR B 


CORR D 


CORR E 


1 


. 1715 


.4936 


.4294 


.4277 


2 


. 1491 


.7761 


. 6292 


. 6328 


3 


. 1767 


. 5375 


.4536 


.4522 


4 


. 1798 


. 6172 


.4354 


.4556 



Standard deviation values x 1 x 10 6 . 



B simulation 
C simulation 
D simulation 
E simulation 



was with fully interactive ocean model, 
was with 10-year mean cycle ocean, 
was with the reduced P . 

was with the improved albedo representation. 
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Several sensitivity runs were made to further explore this 
observation. The range of the albedos was systematically 
increased until values well outside the normal range were 
prescribed. This was done to determine just how insensitive 
the model was to albedo modifications. These runs used a 
single year of integration (1977) in the same manner as for 
the ice strength sensitivity experiment. The surface albedo 
representations examined were as follows: 



Al. 


Albedo code as in 


Semtner 


(1987) . 








A2 . 


New 


albedo code as 


described above 


and 


used in 


Ross 




and 


Walsh (1987) . 














A3. 


All 


ice and 


snow 


albedos 


in 


the 


new 


albedo 


code 




increased by 


0.05. 














A4 . 


All 


ice and 


snow 


albedos 


in 


the 


new 


albedo 


code 



increased by 0.1 and open water albedo increased to 

0.15. 

- A5 . Ice albedo held constant at 0.65, snow albedo held 

constant at 0.80 and open water albedo at 0.15. 

- A6. Ice albedo held constant at 0.85, snow albedo held 

constant at 0.90 and open water albedo at 0.15. 

- A7 . Ice albedo held constant at 0.35, snow albedo held 

constant at 0.40 and open water albedo at 0.10. 

Several other runs were conducted using the high and low 
albedo representations of runs A6 and A7 but with different 
portions of the thermodynamic and dynamic forcing in the 
linked model removed. This was done to determine which of the 
features of the fully linked model reduced its sensitivity to 
albedo changes as compared to previous thermodynamic ice 
models. This work was done in conjunction with the objectives 
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of Chapter VII. For reference, a brief description of the 
various dynamic model variations are noted below. Further 
descriptions of the various runs D1 thru D8 and the rationale 
for their selection are contained in Chapter VII. 



- Dl. 


Base model using the 10-year average ocean. Start 
from Dec 1976 and use 1977 forcing over a single 
integration year. 


- D2. 


Constant ocean heat flux into mixed layer (2.0 x 10’ 
5 deg C/sec, equivalent to approximately 25 W m’ 2 ) . 


- D3 . 


No ocean heat flux. 


- D4 . 


No ocean current stress (level two (40 m) ocean 
velocities = 0.0 m/s). 


in 

Q 


Ocean heat flux = 2.0 x 10' 5 and no ocean current 
stress . 


- D6 . 


No wind stress (surface drag coefficient = 0.0). 


- D7 . 


No wind or ocean current stress. 


- D8 . 


No wind stress, no ocean current stress and no ocean 
heat flux. 



In cases D2 through D8, the description for each run indicates 
the changes from the Base Model (run Dl) . 

B. RESULTS 

The ice thickness distributions for September 1977 from 
runs A6 and A7 are shown in Figures 6.1 and 6.2. These two 
runs had the largest differences in albedos and September was 
the month which showed the greatest differences in ice 
thicknesses. The differences between runs A6 and A7 (and runs 
A1-A5 as well) were small. The ice cover did in general 
increase when the albedo was changed from a low to a high 
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Figure 6.1 



Simulated ice thickness contours (cm) 
September, 1977. High albedo case (run A6) . 
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Figure 6.2 



Simulated ice thickness contours (cm) 
September, 1977. Low albedo case (run A7 ) . 



for 
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value; however, the ice edge advance rarely exceeded 50 km 
seaward. In fact, there were numerous examples in which the 
position of the ice edge and contours within the high 
thickness gradient regions were indistinguishable between the 
high and low albedo cases except for the smallest of 
irregularities . 

Differences between simulations using modified albedos 
became more noticeable in runs D1-D8 . The common differences 
were a southward extension of the marginal ice zone (MIZ) all 
around the Arctic Basin, a southward extension of the EGC, 
increased ice in the Kara Sea and in several cases an increase 
in ice poleward of Spitsbergen. Summer central ice pack 
thicknesses decreased very slightly (approximately 2%) in 
several of the high albedo runs presumably because the thicker 
ice around the margins reduced the amount of ice compression 
in the central pack. Differences between high and low albedo 
simulations were largest in the spring during periods of large 
scale melting. Once the melt was in full progress (and during 
the subseguent freezing period) , the ice fields evolved almost 
identically. 

Figures 6.3 - 6.8 show the ice thickness distributions in 
July 1977 for the high and low albedo conditions from runs D6, 
D7 and D8 . These were the runs in which wind stress was 
reduced to zero. In these cases much more ice was evident in 
the eastern Arctic during summer and fall using the large 
albedo values. Increased southward extension of the MIZ was 
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Figure 6.3 



Simulated ice thickness contours (cm) for July, 
1977. High albedo case (run D6) . 



152 




Figure 6.4 Simulated ice thickness contours (cm) for July, 
1977. Low albedo case (run D6) . 
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Figure 6.5 



Simulated ice thickness contours (cm) for July, 
1977. High albedo case (run D7). 
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Figure 6.6 



Simulated ice thickness contours (cm) for July, 
1977. Low albedo case (run D7 ) . 




Figure 6.7 Simulated ice thickness contours (cm) for July, 
1977. High albedo case (run D8) . 



156 




Figure 6.8 Simulated ice thickness contours (cm) for July, 
1977. Low albedo case (run D8). 
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also observed in the western Arctic but relative to the 
eastern Arctic increase, this change was relatively small. 

In contrast to the small differences in ice cover between 
simulations when the albedo was changed, the differences 
between simulations when the dynamic or thermodynamic forcing 
was changed were dramatic. These effects will be examined in 
Chapter VII. 

C. DISCUSSION 

Prior to this experiment, it was tempting to assume that 
the simple albedo representation used by Semtner (1987) and 
initially used here provided surface albedo values which were 
too low. For a model sensitive to albedo changes, this would 
cause earlier and increased melting and perhaps account for 
the excessive retreat of the ice edge in summer. A higher 
albedo, produced by a more "accurate" albedo representation, 
would reduce the excessive melting and enhance freezing. The 
resulting increased ice thicknesses would provide additional 
resistance to earlier meltoff. 

The results of this experiment showed that this was not 
the case. Increasing the ice and snow albedo did result in 
nominally more ice area and thickness, especially near the 
edges of the pack ice, and less ice was produced at reduced 
albedos. However the differences were small even when the 
albedos were set well outside normal values. 
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One explanation could be that the upper ocean controls the 
extent of the ice cover directly. With a fixed albedo of 
0.10, the surface water absorbs much more heat than the ice. 
Wherever open water and ice are adjacent during spring and 
summer, the warmed water is free to advect under and around 
the ice to promote melting regardless of the albedo of the 
ice. The reduced ice concentration in the MIZ would further 
promote this process. However, this mechanism would also be 
possible in the earlier thermodynamic-only ice models. Since 
those were the models which displayed considerable albedo 
sensitivity, this explanation can be rejected. 

A second explanation is that the inclusion of the 
interactive ocean could permit a feedback mechanism which 
largely compensates for any changes to the frozen-surface 
albedo. For example, the greater ice production rate at 
higher albedos could initiate greater vertical mixing through 
salt extrusion and thereby bring more heat to the surface, 
slowing the ice growth. Similarly the reduced albedo would 
initiate more melting, create a stronger halocline and 
therefore less vertical heat advection. 

A third explanation, and probably the most reasonable, is 
that the ocean heat flux, ocean currents and wind stress are 
much more dominant in controlling the ice edge than the 
atmospheric heat flux changes to the ice resulting from albedo 
modifications. 



159 



The primary difference between this ice model, which has 
been shown to be insensitive to ice surface albedo changes, 
and previous models which were sensitive to albedo changes was 
the inclusion of a fully prognostic ocean. Runs D2 - D8 were 
used to examine the impact of the dynamic features on the 
model's sensitivity to frozen-surface albedo changes. The 
greatest sensitivity was observed in the runs in which the 
wind stress was reduced to zero (D6-D8) . This suggests that 
dynamic forcing, and, in particular the wind stress might act 
to remove, compensate or dominate most of the effect of 
changing the frozen-surface albedo. However the model used 
by Hibler (1980) accounted for the effects of wind stress in 
simulating the annual cycle of Arctic ice cover and still 
found that the model was sensitive to small changes in the ice 
albedo. The major difference remaining between the Hibler 
model and the model used here was that Hibler used a constant 
vertical ocean heat flux. One must therefore assume that it 
is this vertical heat flux from the ocean, which simulations 
from this work indicate varies dramatically in time and space, 
that reduces this model's sensitivity. It follows then that 
the vertical ocean heat flux is likely the most dominant 
factor in determining the thermodynamic balance near the ice 
edge . 

The sensitivity of runs D6 and D7 to albedo changes, 
despite using a simulated variable ocean heat flux, indicated 
that another factor was also important to albedo sensitivity. 
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Ice thicknesses and areal coverage were severely reduced in 
run D6 compared to run D7, and run D6 also displayed more 
sensitivity to the albedo change. The areas of thin ice 
became ice free much more quickly and did not extend as far 
when the albedo was low. It would therefore appear that the 
model's albedo insensitivity was also dependent on the ice 
thickness . 

Further support for this is evident by comparing runs D3 
and D8 . In both cases the ocean heat flux was reduced to 
zero; however D3 retained the dynamic forcing from wind and 
ocean currents. The lack of dynamic forcing in D8 reduced the 
simulated ice thicknesses in the central pack and the 
thickness gradients in the MIZ. The extended, thin ice MIZ 
in D8 was much more susceptible to the albedo change than the 
strong thickness gradient MIZ in D3 . 

There is an apparent contradiction in that D7 has thicker 
ice than D6 despite reduced dynamic forcing. However it 
appears that the water stress, which is accounted for in Run 
D6 but not in D7 , advects the ice away from the thick 
compression areas to regions of high ocean heat flux where it 
is melted. This reduces the total ice volume and ice 
thicknesses, especially in the central pack and along the 
Canadian boundary of the Beaufort Gyre. This process may also 
explain the reduced ice thicknesses immediately adjacent to 
the Canadian Archipelago and northern Greenland evident in the 
ice strength experiment (Chapter V) . 
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D. CONCLUSIONS 



The use of only one integration year vice several years or 
decades of integration obviously limited the degree of ice 
cover modification possible between model runs with different 
surface albedo representations. However the objective here 
was not to determine what final end condition of ice cover 
could be produced but whether or not alteration of the albedos 
within reasonable limits could improve the simulation of 
interannually varying ice cover. 

The largest error in the simulated ice fields was the 
excessive summer melt. A secondary, but still important 
error, was excessive winter ice coverage in the Barents Sea. 
This model appears quite insensitive to albedo changes of the 
snow/ice cover. Surface albedos had to be set well above 
reasonable levels to produce even a nominal decrease in the 
meltback in the western Arctic and the same albedo changes 
marginally worsened the simulation of the ice edge in the 
Barents Sea. Therefore, changing the albedo representation 
was not an effective method of improving the accuracy of the 
simulated ice cover. In fact, the insensitivity of this model 
to albedo modifications would suggest that even the very 
simple albedo representation initially used would be 
sufficient and would not degrade the model's accuracy. 

The primary reason for this model's insensitivity to 
albedo changes in the frozen surface appears to be the 
dominance of the vertical ocean heat flux in the thermodynamic 
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balance near the ice edge. Albedo sensitivity in the model 
is only apparent in those few, small regions where both the 
vertical ocean heat flux is small and the ice is thin. Under 
these conditions, the changes to the thermodynamic balance 
induced by albedo changes in the ice cover are large enough 
to cause noticeable changes in the ice concentration and 
extent. 

It should be reiterated here that this study deals with 
the ice cover on an average monthly basis at a fairly large 
scale. Albedo does have a major impact on the ice cover on 
shorter time scales as evidenced, for example, by the flash 
melt phenomenon. The entire ice surface in a region can melt 
in a few days when the reduced albedo of the melting surface 
accelerates the melting process. However, on the larger time 
and space scales used here, this type of effect is largely 
smoothed out leaving the consistent longer term ocean heat 
flux as the dominant controller of the ice cover. 

The ocean currents, although not a major factor in albedo 
sensitivity, play a significant role in the evolution of ice 
thickness. They appear to advect the ice away from the thick 
compression regions to regions of high heat flux where the ice 
is usually melted. This has the effect of reducing the ice 
thickness in the compression regions, particularly along the 
northern limit of the Canadian Archipelago. 
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VII 



DOMINANT MECHANISM EXPERIMENT 



The albedo sensitivity experiment included several runs 
in which various dynamic and thermodynamic processes were 
removed from the base linked model to determine which 
processes were primarily responsible for the insensitivity of 
the model to albedo changes. Although the albedo changes 
continued to have relatively little impact on the simulated 
ice fields, the total ice field changed considerably from the 
base model output as the different processes were omitted. 

This experiment was conducted to determine which elements 
of the thermodynamic and dynamic forcing were dominant in the 
linked ice-ocean model. 

A. EXPERIMENTAL SETUP 

The same base model as used in the albedo sensitivity 
study was used for this experiment. A single integration year 
was considered sufficient due to the high sensitivity of this 
model to dynamic changes. The year 1977 was again chosen as 
the integration year because the previous work had provided 
a large base of simulated 1977 ice field data for comparison. 
The following sensitivity runs were examined: 

- Dl. Base model. 10-year average ocean. Starting from 

Dec 1976 output fields. 1977 forcing. Single 
integration year. 

- D2 . Constant ocean heat flux into mixed layer (2.0 x 

10’ 5 deg C/sec, equivalent to approximately 2 5 W/m 2 ) . 
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- D3 . No ocean heat flux. 

- D4. No ocean current stress (level two (40 m) ocean 

velocities = 0.0 xn/s) . 

- D5 . Ocean heat flux = 2.0 x 10’ 5 and no ocean current 

stress. 

- D6. No wind stress (surface drag coefficient = 0.0). 

- D7. No wind or ocean current stress. 

- D8 . No wind stress, no ocean current stress and no ocean 

heat flux. 

In cases D2 through D8, the description for each run indicates 
the changes from the Base Model (run Dl) . 

B. RESULTS 

The winter (April) and summer (August) ice thickness 
contours from run Dl are shown in Figures 7.1 - 7.2 for 

comparison with the contours from runs D2-D8. 

The ocean heat flux value chosen for run D2 was the 
approximate median of the simulated average annual flux cycle 
over the Arctic region. Applying this constant value of 25 
W/m 2 over the entire grid heavily exaggerated the basin-wide 
simulated ice cycle. The winter ice cover expanded much too 
far south into the Norwegian Sea and the summer ice cover was 
reduced to a small fraction of the observed cover (Figures 
7.3 - 7.4). The central ice pack normally receives little or 
no ocean heat flux. The relatively large heat flux applied 
in this run rapidly reduced the pack ice thickness and areal 
coverage. By August, the central pack was almost completely 
melted off. The subsequent freezing in early winter produced 
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Figure 7.1 Simulated ice thickness contours (cm) for April, 
1977. Base model (run Dl) . 
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Figure 7.2 Simulated ice thickness contours (cm) for August, 
1977. Base model (run Dl) . 
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Figure 7.3 



Simulated ice thickness contours 
1977. Constant ocean heat flux 
D2 ) . 



(cm) for April, 
=25 w/m 2 (run 
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Figure 7.4 



Simulated ice thickness contours (cm) for August, 
1977. Constant ocean heat flux = 25 W/m 2 (run 
D2 ) . 
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a totally unrealistic simulated winter ice cover. Although 
the prescribed heat flux in the Barents and Norwegian Seas 
during winter was positive, it was still considerably less 
than in the base model case (run Dl) . This resulted in 
considerably less melting and extended ice coverage there. 

Run D3 used an ocean heat flux set to zero. This produced 
even heavier ice during winter in the North Atlantic and 
Greenland Sea (Figure 7.5) because the ocean heat flux in that 
region was reduced even further from the values used in run 
Dl. The ice also remained longer and stayed thicker around 
the western ice boundary where excessive meltback occurred in 
run Dl in summer (Figure 7.6) . However the meltback was still 
somewhat severe compared to the observed data. For the 
purposes of these observations the western ice boundary 
included the Beaufort, Chukchi and East Siberian Seas. The 
central pack thickness and distribution remained relatively 
unchanged from the base model as the normal simulated heat 
flux there was also zero. 

Run D4 had no stress from ocean surface currents. In this 
case the pack ice thicknesses decreased slightly but the 
relative distribution of the ice thickness changed very 
little. In comparison to the base run, the ice cover over the 
Beaufort Gyre was rotated counter-clockwise approximately 
fifteen degrees (Figures 7.7 - 7.8). The MIZ expanded outward 
in all regions except the EGC where the ice advection was 
suppressed. The thickness gradient marking the poleward side 
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Figure 7.5 Simulated ice thickness contours (cm) for April, 
1977. Ocean heat flux = 0.0 W/m 2 (run D3). 
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Figure 7.6 Simulated ice thickness contours (cm) for August, 
1977. Ocean heat flux = 0.0 W/m 2 (run D3 ) . 
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Figure 7.7 Simulated ice thickness contours (cm) for April, 
1977. Level two current velocities =0.0 m/s (run 
D4 ) . . 
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Figure 7.8 Simulated ice thickness contours (cm) for August, 
1977. Level two current velocities = 0.0 m/s (run 
D4 ) . 
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of the Transpolar Drift Stream (TPD) also became less intense. 
Three supplementary runs were conducted with the same dynamic 
conditions as run D4 but using different ice strength values 
as in the ice strength experiment. As ice strength was 
increased, which reduced the plasticity of the ice, the main 
features separating run D4 from run D1 became less noticeable. 

In run D5, a similar counter-clockwise rotation to the ice 
field as in run D4 was observed for the Beaufort Gyre ice 
cover. The reduced advection in the TPD and EGC was evidenced 
by the reduced thickness gradients in the current regions (See 
Figure 7.9). The combined reductions in ice pack thickness 
from higher ocean heat flux in the central Arctic and no 
dynamic ocean forcing resulted in a complete ice meltback by 
August . 

The direct dynamic forcing of the wind on the frozen 
surface was removed in run D6 by setting the ice/air and 
snow/air drag coefficients to zero. Ice thickness was reduced 
by approximately 50% within six months and the thickness 
distribution pattern also changed substantially (Figures 7.10- 
7.11). The central pack shifted over 1000 km towards the 
Beaufort Sea, but due to the reduced ice thickness, the 
thickness gradient there changed very little. Ice remained 
in the western boundary seas at least one month longer than 
the standard case and the initial MIZ retreat in that region 
was not as severe. However the overall reduced thickness 
resulted in a smaller total ice cover in the Arctic by August. 
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Figure 7.9 Simulated ice thickness contours (cm) for April, 
1977. Ocean heat flux and level two velocities 
= 0.0. (run D5) . 
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Figure 7.10 Simulated ice thickness contours (cm) for April, 
1977. Wind stress = 0.0. (run D6) . 
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Figure 7.11 Simulated ice thickness contours (cm) for August, 
1977. Wind stress = 0.0. (run D6) . 
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The advection of ice in the TPD and EGC was noticeably 
reduced. 

The run D7 simulation removed all direct dynamic forcing 
from both the air and ocean. No horizontal ice velocity was 
permitted. The same large shift of the central ice pack 
towards the Beaufort Sea was observed similar to run D6 but 
the ice thicknesses were not as severely reduced. The MIZ 
retreat in the Beaufort Sea started earlier but did not move 
as far poleward as in run D6. Run D7 showed a noticeable ice 
retreat poleward of Spitsbergen and in the EGC which would 
correspond to the removal of the TPD and EGC (Figures 7.12 - 
7.13) . A comparison of the summer thickness contours for runs 
D5, D6 and D7 shows that in those simulations where the 
central pack is reduced to thicknesses on the order of 150 cm 
or less, the subsequent evolution of the ice cover changes 
dramatically. There appears to be a critical thickness below 
which the influence of the various mechanisms controlling the 
ice changes significantly. 

Finally run D8 used the same dynamic limitations as run 
D7 ; however the ocean heat flux was also set to zero. The 
same excessive ice cover over the Barents Sea as was observed 
in run D3 was repeated and the same general ice thickness and 
distribution pattern as for run D7 was produced. However in 
contrast to all the other runs, the ice remained over the 
western boundary seas for the entire year-long simulation. 
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Figure 7.12 Simulated ice thickness contours (cm) for April, 
1977. Wind and ocean current stress = 0.0. (run 
D7 ) . 
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Figure 7.13 Simulated ice thickness contours (cm) for August, 
1977. Wind and ocean current stress = 0.0. (run 
D7 ) . 
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Furthermore, the ice edge did not retreat poleward of 
Spitsbergen, a condition closer to the observed fields than 
the previous simulations (Figures 7.14 - 7.15). 

C. DISCUSSION 

This experiment was conducted primarily to determine 
whether ocean heat flux or water stress from ocean currents 
was the more important process determining the limit and 
thickness distribution of the ice cover. These are the two 
primary processes which become more realistic once the 
prognostic ocean is included in the linked model. The results 
indicated that such a simple determination was not possible. 
As has been noted in the previous chapters, there appear to 
be many factors working in combination and often with opposing 
effects. 

The direct dynamic influence of ocean surface currents 
influences the simulated ice cover by distorting it in the 
direction of the surface currents. The Canadian Basin ice 
pack is rotated clockwise by the Beaufort Gyre and ice is 
advected along the streamlines of the TPD and the EGC. This 
increases the thickness gradient poleward of the TPD and 
pushes the summer ice edge in the eastern Arctic towards 
Spitsbergen and well south along the east coast of Greenland. 
The region of thickest multi-year ice is shifted towards the 
Fram Strait. The ocean currents converge the ice to the 
center of the gyre and into a region off the north coast of 



182 




Figure 7.14 Simulated ice thickness contours (cm) for April, 
1977. Wind and ocean current stress and ocean 
heat flux = 0.0. (run D8) . 
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Figure 7.15 Simulated ice thickness contours (cm) for August, 
1977. Wind and ocean current stress and ocean 
heat flux = 0.0. (run D8). 
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Greenland and the Canadian Archipelago. This accounts for at 
least a portion of the thicker ice evident in those regions. 

The ice is advected parallel to the coastline all around 
the perimeter of the Canadian Basin. This reduces the ice 
thickness immediately poleward of the Canadian Archipelago 
but the ice does not collect in the Beaufort Sea since most 
of it is quickly melted by the high ocean heat flux there. 
The ice that is left continues to be pushed along the 
coastline out of the Beaufort and Chukchi Seas and into the 
Laptev Sea. 

The ice cover in the Norwegian and Barents Seas is nearly 
unaffected when the water stress from ocean currents is 
ignored. However changes to the ocean heat flux result in 
major changes to the simulated ice cover. Of course the heat 
flux at any position is a function of earlier horizontal and 
vertical advection. Therefore the currents do serve a purpose 
in this region. However it is notable that the dynamic impact 
of these currents at the space and time scales used here is 
small . 

The wind stress appears to have a much stronger effect on 
the ice cover than the ocean current stress at these spatial 
and temporal scales. The ice thickness and distribution is 
drastically altered when wind stress is removed from the 
simulations. The wind stress appears to provide the majority 
of the compressive forcing which produces the realistic ice 
thicknesses in the central pack and convergence regions. For 
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the 1977 forcing year, the wind provided an offshore component 
to the ice edge position around the entire western half of the 
Arctic Basin, particularly in the East Siberian Sea. In the 
eastern Arctic Basin the wind stress appeared to increase the 
ice advection along the TPD and the EGC and the ice cover was 
pushed southward in the Barents Sea. It should be noted that 
the Walsh monthly mean wind fields used in this study vary 
interannually in direction and intensity (Walsh, private 
communication) . Therefore considerable differences in the 
wind forcing and corresponding changes to the ice fields would 
probably be observed had a different forcing year been used 
in this experiment. Nevertheless, the relative strength of 
this forcing mechanism would be expected to remain similar. 

The ocean heat flux, wind stress and ocean currents form 
a feedback which amplifies the influence of each of the 
individual mechanisms on the melt cycle. For example, the 
ocean heat flux thins the ice which then allows the dynamic 
forcing to have more effect. Conversely, the thicker ice 
regions which develop over areas of reduced or negative ocean 
heat flux become less sensitive to the dynamic forcing. This 
feedback is stronger when the ice strength is reduced as 
described in the ice strength experiment. 

D. CONCLUSIONS 

The simulations conducted in this experiment support the 
view that the importance of the various dynamic and 
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thermodynamic mechanisms controlling the monthly variations 
of ice coverage vary regionally across the Arctic Basin. The 
ocean heat flux appeared to be the dominant mechanism 
controlling the extent of the monthly average ice cover while 
the wind and current stresses appeared to be the mechanisms 
determining the ice thickness and distribution. However the 
dynamic forcing from winds and currents had quite different 
effects in the different regions of the Arctic Basin. In the 
western basin which was incorporated in regions 1 and 2, the 
wind and currents acted in concert with the ocean heat flux 
to clear the ice from the boundary seas. In contrast, these 
same forces acted to increase the ice cover in the TPD and EGC 
(region 3) while in the Barents and Norwegian Sea (region 4) , 
the wind stress tended to oppose the effects of ocean heat 
flux and current stress had little impact. In general, the 
forcing from ocean currents was usually weaker than from ocean 
heat flux or wind stress. The exceptions were in the strong 
current regions along the coast in the Canadian Basin and 
along the axis of the TPD and EGC. Inclusion of wind and 
current stresses always tended to increase the ice area, 
either directly or indirectly. In regions 3 and 4 the effect 
was direct by off-ice advection. In regions 1 and 2 the 
effect was indirect. The dynamics initially caused the ice 
to converge which thinned the ice on the boundaries thereby 
allowing more ice growth and subsequent convergence. As the 
central pack became thicker, it also became stiffer. The 
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thickness gradients began to extend further from the central 
core and the MIZ also increased in thickness. The ocean flux 
was unable to melt the same extent of a thicker MIZ; therefore 
the ice area expanded. 

The ocean does not appear to control the ice coverage 
significantly through direct dynamic forcing. However it does 
appear to have a very important thermodynamic role. The three 
dimensional ocean circulation positions the heat contained in 
the intermediate depth ocean layers in the appropriate 
regions. This heat is then available whenever the vertical 
circulation, driven by conditions at the surface, extends deep 
enough to tap it. The ocean essentially preconditions the 
ocean below the mixed layer so that the heat necessary to 
control the ice edge at any location is available. 

The dramatic changes to the ice cover evident in the 
simulations where wind stress and ocean heat flux were 
modified emphasized the sensitivity of the model to these 
parameters. The importance of realistic ocean heat flux was 
already suggested in previous chapters and is further 
supported by the results of this experiment. This experiment 
also emphasized the importance of appropriate wind forcing 
fields if the monthly average ice cover is to be simulated 
accurately in the Arctic. 
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VIII. STATISTICAL CONSIDERATIONS 



The time series of "observed" ice concentration data 
covered the same regions at the same grid scale and over the 
same time period as the simulated data. This permitted 
calculation of the correlations between the various model 
simulations and the observed fields. The differences between 
the correlations calculated for each new experiment were 
compared to determine if the changes incorporated into the 
model improved the accuracy of the simulated ice concentration 
fields. A large increase in the correlation coefficient 
indicated that the modifications made to the model improved 
the accuracy of the output. The sign of the correlation 
differences were usually consistent indicating that the 
modifications made had produced a noticeable and consistent 
change. However statistical analysis of the results was 
required in order to determine if these changes were 
statistically significant. 

Chervin and Schneider (1976) and Chervin (1980) 
investigated a similar problem as applied to atmospheric 
general circulation models (GCMs) . They suggested that a 
knowledge of the "noise climatology" or natural variability 
inherent in a GCM was a prerequisite for judging the 
significance of the results from prescribed change model 
experiments. Change model experiments are similar to what has 
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been done in this work in which various elements of a base 
model are changed to investigate changes in the simulated end 
conditions. Chervin and Schneider (1976) noted that the basic 
problem was one of 

distinguishing between signal (that part of the difference 
between the results of a prescribed change experiment and 
an unperturbed control case which is attributable to the 
prescribed change) and noise (some measure of inherent 
model variability) in a prescribed change response. 

The model variability, in the case of the atmospheric 
GCM's referred to above, was a result of the numerical 
representation of random processes in the atmosphere. Each 
run of such a model allows for some randomness and as a result 
the simulated fields will vary from run to run with no change 
in boundary conditions or numerical code within the model. 
The numerical model used here does not have a similar output 
variability. The simulated ice fields produced will not 
change unless the boundary conditions or numerical code is 
modified. However in this work, it is not the model output 
fields which are directly compared but the difference between 
models of the correlations between simulated and observed ice 
concentration fields. Therefore it is these correlation 
differences which must be examined to determine their 
statistical significance. 

A. METHOD 

The correlations calculated in this work use 120 sample 
pairs (every month for ten years) . Each value in each paired 



190 



sample was the average of several hundred grid values of ice 
area contained within each of the four designated regions. 
A common method to establish the variability and subsequently 
the confidence intervals of the correlation differences, 
required that subsets of the paired samples be selected and 
assumed independent. The time series of 120 paired samples 
were divided into subsets which were deemed independent and 
correlations calculated for each. Means and standard 
deviations for these subset correlations were then determined. 
However if the number of subsets was too small or if the 
subsets were not independent of each other, their means were 
likely to be biased. Significant persistence of ice 
concentration anomalies has been observed in areas of the 
eastern Arctic to lags of several months and in exceptional 
cases to over a year (Fleming, 1987) . The method presented 
here strikes a compromise balancing the opposing requirements 
of independence versus a fairly large number of subsets. A 
subset size of one year was selected because it was 
exceptional to find dependence at lags over a year. This 
determined that the maximum number of subsets was 10. 

The standard formulation used to calculate the 95% 
confidence interval was: 

XB ± t 0-975 x SD x (n) 1/2 

where XB was the average correlation over 10 subsets, t was 
the Student t distribution value at 95% confidence level and 
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9 degrees of freedom, SD was the estimated standard deviation 
of the 10 subset correlations and n was the number of years 
( 10 ). This formula was derived from normal theory which 
required n > 20 and observations independent and identically 
distributed. Obviously our sample size of 10 made this 
technique suspect, however it was the largest sample possible 
in view of the dependence consideration. 

An alternative technique used to investigate the 
statistical significance of the difference between the various 
correlations was called "Jackknife." (For a complete 
description see, for example, Mosteller and Tukey (1977)). 
This technique offers a way to establish sensible confidence 
limits in complex situations in which the observations are not 
necessarily independent and identically distributed. The 
correlations were calculated with all the data and then after 
dividing the data into subsets, correlations were recalculated 
by leaving out each subset, one at a time. The subset size 
chosen was again one year. Pseudo values were determined by 
weighting the correlations appropriately to account for the 
sampling "overlap" and the Student t tables were then used to 
establish the confidence limits. The procedure is outlined 
below. 

The correlations calculated for two different models using 
all 10 years of data are r,(0) and r 2 (0). The correlation 
difference, Y(0), is then 

Y ( 0 ) = ri (0) - r 2 ( 0) 
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If r,(n) and r 2 (n) are the correlations leaving out year n, 
where n=l , 2 , . . . , 10 , then 

Y (n) = r,(n) - r 2 (n) 

The n th pseudo value, Y (n) , is calculated as follows: 

Y* (n) = (10 x Y (0) ) - (9 x Y (n) ) 

where 10 and 9 are the weights appropriate to our one year 
subset size. The variances, s 2 , of the ten pseudo values are 
calculated and the approximate two-sided 95% confidence limit 
for the difference of the true correlation was given by: 

Y ( 0) ± t 0 975 x ( s 2 /10) 1/2 

where t 0 975 = 2.2 6 2. This value was found in the Student t 
table using 9 degrees of freedom. 

B. RESULTS 

The common and Jackknife techniques were applied to the 
correlations calculated in the Interactive Ocean Experiment 
(Table 5) . For the common method, means, standard deviations 
and confidence intervals were calculated from each group of 
ten subset correlations. These are shown in Table 8. 

Confidence intervals were calculated for both sets of 
correlations (the B series from the model with the interactive 
ocean and the C series from the model using the ten year 
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TABLE 8 



MEANS, STD DEVIATIONS AND CONFIDENCE INTERVALS OF THE 
CORRELATIONS OF 10, 1-YEAR SUBSETS. RESULTS FOR 
EXPERIMENT B (FIRST LINE) AND C (SECOND LINE) 

IN ALL FOUR REGIONS 



ICE AREA TIME SERIES 



REGION 


MEANS 

CORR B/CORR C 


STD DEV 
CORR B/CORR C 


95% CONF. INT. 
CORR B/CORR C 


1 


.9322 


.0580 


.8907-. 9737 




.9200 


.0504 


.8839-. 9561 


2 


.9323 


.0529 


.8945-. 9701 




.9198 


.0627 


.8750-. 9646 


3 


.9628 


.0164 


.9511-. 9745* 




. 9396 


.0271 


.9202-. 9590* 


4 


.9500 


.0190 


.9364-. 9636* 




.9061 


.0466 


.8728-. 9394* 



ANOMALY TIME SERIES 



REGION 


MEANS 

CORR B/CORR C 


STD DEV 
CORR B/CORR C 


95% CONF. INT. 
CORR B/CORR C 


1 


. 3472 


.4404 


.0322-. 6622 




.2779 


.3864 


. 0015-. 5543 


2 


.7198 


. 1696 


.5985-. 8411* 




. 5670 


.3614 


.3085-. 8255 


3 


.4021 


.2422 


.2289-. 5753 




.3305 


.2583 


. 1457- . 5153 


4 


. 6266 


.2428 


.4529-. 8002 




.4952 


.2354 


.3268-. 6636 



B simulation was with fully interactive ocean model 
C simulation was with 10-year mean cycle ocean 
* difference between the correlations is 95% significant 
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average ocean condition) . The intervals were then compared 
to determine if they confirmed each other and how important 
the variance of the correlations was. 

The differences between the correlations for regions 3 and 
4 of the ice area time series were determined to be 
significant at the 95% confidence level. The average of the 
C correlations fell outside the confidence limits of the B 
correlations and the average of the B correlations fell 
outside the confidence interval of the C correlations. This 
provided some measure of confirmation. The correlation 
differences between B and C for the ice area anomaly time 
series were primarily insignificant at the 95% confidence 
level. One exception was the B model's region 2. In that 
case, the C model average correlation fell outside the B 
model's correlation confidence interval. However, the 
opposite was not true due to a much higher correlation 
variance for the C model. 

In the process of applying the Jackknife technique, the 
correlation differences (B subtract C) for all four regions, 
for both time series, and for 10 pseudo value subsets were 
calculated. These are shown in Table 9. Note that they were 
consistently positive. Table 10 shows the calculated pseudo 
values, their means, their standard deviations and the 
calculated 95% confidence intervals for the ice area time 
series. Table 11 shows the same information for the ice area 



195 



TABLE 9 



CORRELATION DIFFERENCES BETWEEN B AND C FOR THE TOTAL TIME 
SERIES AND FOR THE 10 PSEUDO SUBSETS (ONE YEAR REMOVED) 



ICE AREA TIME SERIES 



PSEUDO 

SUBSET 


REGION 1 


REGION 2 


REGION 3 


REGION 4 


TOTAL 


. 0274 


. 0464 


.0355 


.0854 


1 


. 0295 


. 0328 


.0388 


. 1079 


2 


. 0266 


. 0453 


.0290 


.0883 


3 


.0340 


. 0494 


. 0388 


.0540 


4 


. 0221 


. 0461 


.0341 


. 0941 


5 


. 0316 


.0579 


.0353 


. 0809 


6 


. 0270 


. 0469 


.0436 


. 0888 


7 


. 0306 


.0533 


. 0362 


.0928 


8 


. 0248 


. 0482 


. 0306 


.0792 


9 


. 0233 


. 0306 


.0322 


. 0825 


10 


. 0233 


. 0494 


.0349 


. 0814 



ICE AREA ANOMALY TIME SERIES 



PSEUDO 

SUBSET 


REGION 1 


REGION 2 


REGION 3 


REGION 4 


TOTAL 


. 0973 


. 1780 


.2336 


.2840 


1 


. 1093 


. 1038 


.2496 


.3731 


2 


. 0922 


. 2539 


. 1826 


.3110 


3 


. 1271 


. 1647 


.2911 


.1577 


4 


. 0974 


. 1978 


.2185 


.3166 


5 


. 1407 


.2125 


.2104 


.2498 


6 


. 0976 


.2034 


. 2978 


.2807 


7 


. 1257 


. 1956 


.2644 


.3094 


8 


. 0819 


. 1714 


. 1586 


.2004 


9 


. 0792 


. 1068 


.2240 


.3455 


10 


.0312 


. 1797 


.2263 


.2971 
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TABLE 10 



PSEUDO VALUES FOR THE ICE AREA TIME SERIES 
MEANS, STD DEVIATIONS AND 95% CONFIDENCE INTERVALS OF THE 

PSEUDO VALUES 



ICE AREA TIME SERIES 



YEAR 

REMOVED 


REGION 1 


REGION 2 


REGION 3 


REGION 4 


1 


. 0085 


. 1688 


. 0058 


-.1171 


2 


. 0346 


.0563 


.0940 


. 0593 


3 


-.0320 


. 0194 


. 0058 


. 3680 


4 


. 0751 


. 0491 


. 0481 


. 0071 


5 


-.0104 


-.0571 


. 0373 


. 1259 


6 


. 0310 


. 0419 


-.0374 


. 0548 


7 


-.0014 


-.0157 


. 0292 


. 0188 


8 


. 0508 


. 0302 


. 0796 


. 1412 


9 


. 0643 


. 1886 


. 0652 


.1115 


10 


. 0643 


. 0194 


.0409 


. 1214 


STD DEV 


.0361 


.0756 


.0388 


. 1246 


MEANS 


. 0285 


. 0501 


.0369 


. 0891 


CONFIDENCE 


.0016 TO 


-.0077 TO 


.0077 TO 


-.0037 TO 


INTERVAL 


.0532 * 


. 1005 


.0633 * 


. 1745 



* indicates that the correlation difference is 95% significant 
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TABLE 11 



PSEUDO VALUES FOR THE ICE AREA ANOMALY TIME SERIES 
MEANS, STD DEVIATIONS AND 95% CONFIDENCE INTERVALS OF THE 

PSEUDO VALUES 



ICE AREA ANOMALY TIME SERIES 



YEAR 

REMOVED 


REGION 1 


REGION 2 


REGION 3 


REGION 4 


1 


-.0107 


.8458 


.0896 


-.5179 


2 


. 1432 


-.5051 


. 6926 


.0410 


3 


-.1709 


.2977 


-.2839 


1 .4207 


4 


. 0964 


-.0002 


.3695 


-.0094 


5 


-.2933 


-.1325 


.4424 


. 5918 


6 


. 0946 


-.0506 


-.3442 


. 3137 


7 


-.1583 


. 0196 


-.0436 


.0554 


8 


. 2359 


.2374 


.9086 


1.0364 


9 


.2602 


.8188 


. 3200 


-.2695 


10 


. 6922 


. 1627 


. 2993 


. 1661 


STD DEV 


. 2794 


.4145 


. 3997 


. 5883 


MEANS 


. 0889 


. 1694 


. 2450 


. 2828 


CONFIDENCE 


-.1026 TO 


-.1185 TO 


-.0523 TO 


-.1368 TO 


INTERVAL 


. 2972 


. 4745 


. 5195 


.7048 
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anomaly time series. The correlations are considered 
significantly different using the Jackknife method if the 
confidence interval does not include zero. This is the case 
for region 1 and 3 of the ice area time series; however, no 
similar cases exist for the anomaly time series. 

C. DISCUSSION 

The consistently positive differences between the 
correlations from the B and C models would suggest that the 
B model was a significant improvement over the C model. 
However neither of the two statistical approaches used here 
was able to confirm this at the 95% confidence level for more 
than three of the eight correlation difference cases. The 
main reason for this was the high degree of correlation 
variability . 

Correlation variability was most noticeable in the ice 
area anomaly time series because the dominating influence of 
the annual cycle was removed, thereby reducing the correlation 
coefficient values and allowing more correlation variability. 
This was also apparent by comparing the actual time series 
plots (Figures 4.19 - 4.22 with 4.31 - 4.38). The variability 
indicates that the model did very well for some years in some 
regions yet very poorly in others. The simulated ice 
condition at any time was a direct function of the forcing and 
initial conditions. The model had no inherent variability. 
Therefore the changing performance of the model had to be a 
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result of some other cause. One explanation could be that 
some mechanism which was not properly represented in the model 
was controlling the ice and varying dramatically in importance 
from year to year, and region to region. The model would 
perform fine when this mechanism had little influence but 
would be unable to account for those periods or regions when 
the mechanism was important. Alternatively, the prescribed 
forcing may be in error. In particular the monthly averaged 
atmospheric forcing may have too long a temporal resolution 
to adequately represent strong, short term episodic changes. 
The ice dynamics experiment indicated that these could be 
important because wind stress was able to alter the ice field 
dramatically in short periods of time and the subsequent 
evolution of the ice field was also modified. Of course some 
combination of model limitation and forcing error was also 
possible . 

D. CONCLUSIONS 

This work was done to determine if the correlation 
differences between the different model simulations and 
observed fields of ice concentration could be shown to be 
statistically significant. In general, it appears that the 
correlations determined between modelled and observed fields 
of ice area are apparently not yet stable enough to do so. 
The improvement in the correlations from model C to model B 
was consistent in all regions, for all subsets and for both 
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the ice area and anomalous ice area time series. This would 



certainly suggest that the improvement was real. However the 
correlations were so variable, particularly in the ice area 
anomaly case, that supporting statistics could not be 
provided . 
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IX. SUMMARY 



This work has involved the development of a linked ice- 
ocean numerical model capable of simulating the annual cycle 
and interannual variations of ice cover in the Arctic. The 
accuracy of the simulated seasonal ice area and ice edge 
position appears to be an improvement over any other large 
scale models operating at present, yet the computational 
requirements have been maintained at very reasonable levels. 

A. MAIN CONCLUSIONS 

The objectives for this thesis have been met. It has been 
shown that the inclusion of an interactive, prognostic ocean 
component in the ice model provides substantial improvement 
to simulations of both the annual cycle and interannual 
variations of the ice cover. Closer examination of the 
effects of the various ice control mechanisms indicated that 
the main reason for improved ice cover simulation with the 
prognostic ocean was a more realistic representation of the 
ocean heat flux. The ocean heat flux varied considerably 
between regions and between seasons. The surface currents 
were also highly variable. However, in general, the forcing 
from ocean currents has a second order effect on the ice 
cover. Exceptions to this were evident in the strongest 
current regions (East Greenland Current, Transpolar Drift 
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Stream and boundaries of the Beaufort Gyre) where the dynamic 
impact of the currents had a more pronounced effect on the 
ice. 

The linked model was sensitive to changes in the ice 
rheology. Based on comparison with highly averaged observed 
thickness data (Bourke and Garret, 1987) the optimum value 
chosen for the strength parameter in the "bulk viscous" 
rheology was P* = h x 10 4 . This value reduced the rigidity of 
the ice pack, allowing greater compression and therefore more 
realistic ice thickness and distribution. However, the 
reduction in the ice strength degraded the correlation between 
simulated and observed total ice cover. The ice retreat in 
the western Arctic, which was already somewhat severe at the 
higher ice strength, receded further when P* was reduced. 
Since there was no compensating improvement in the winter ice 
extent, the correlation coefficients decreased. 

The increased melting in the western Arctic at the 
reduced ice strength appeared to be a result of a feedback 
between ice formation rates and ocean heat flux. The reduced 
ice strength permitted a stronger response of the ice to 
dynamic forcing. This forcing was generally off-shore in the 
western Arctic. The ice was forced to converge in the central 
pack and along the north shore of the Canadian archipelago, 
opening up large areas of open water and thin ice over the 
shelves in the western Arctic seas. This permitted faster 
cooling of the surface water, greater ice production, more 
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salt extrusion and greater vertical mixing. The mixing was 
strong enough to break through the halocline and increase the 
ocean heat flux to the surface. This resulted in more 
extensive melting. 

The fully linked model with the prognostic ocean was quite 
insensitive to albedo changes in the ice cover. This result 
was a clear contrast to previous model studies which 
demonstrated strong sensitivity to this parameter. The main 
mechanism included in this model but not in the models which 
displayed albedo sensitivity was the interannually variable 
ocean heat flux. This term appears to dominate the 
thermodynamic balance near the ice edge. The ocean heat flux 
thereby controls the concentration of the ice in this region 
and the overall extent of the ice pack. The thermodynamic 
impact of changing surface albedo is relegated to minor 
importance . 

Results from the dynamic mechanism experiment supported 
the idea that the ocean heat flux dominated the albedo effect. 
However these results also indicated that there were 
conditions, particularly when the ice was thin, when the 
albedo did become important. Those cases which showed some 
albedo sensitivity had a combination of a small ocean heat 
flux and thin ice. Under these conditions, the change in the 
thermodynamic balance due to a change in the surface albedo 
could become relatively large and the model's albedo 
sensitivity was increased. These conditions were possible 
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over the shallow shelves, in newly frozen leads within the 
central pack and along some areas of the MIZ. Realistic 
representation of the frozen-surface albedo was therefore 
considered worthwhile and the improved albedo scheme of Ross 
and Walsh (1987) was incorporated into the model. 

Examination of the relative importance of the various 
dynamic and thermodynamic ice cover forcing mechanisms 
indicated that each part of the fully linked model had a 
unique but interactively important role in the evolution of 
the ice cover. The ocean heat flux appeared to be the overall 
dominant factor controlling the position of the ice edge and 
the extent of the ice cover. Within the polar pack, it was 
the dynamic forcing, and in particular, the wind forcing which 
controlled the ice thickness and thickness distribution. 

The dynamic portion of the ice model can be viewed as 
working in concert with the ocean model to produce the desired 
ice cover. The ocean model does not provide much of a dynamic 
effect from surface currents but it does provide an important 
thermodynamic control. The ocean circulation below the mixed 
layer acts to position heat below those regions where 
observations indicate that the ice is usually thin or melted 
away. The dynamic ice model in turn tends to compress the ice 
in those areas where it should regionally be thickest and thin 
the ice in those areas where the heat flux will be used to 
melt back the ice cover. The linkage between these two 
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processes is controlled by conditions at the surface and the 
response of the mixed layer. 

Vertical mixing needs to be induced to tap the 
intermediate layer heat source while stratification must be 
induced to shut the heat off. However, the large horizontal 
diffusion coefficient required for numerical stability 
probably causes excessive lateral mixing of heat within the 
intermediate layer. The increased heat available in some 
regions may be sufficient to maintain excessive heat flux into 
the mixed layer, despite stratification from fresh meltwater. 

It must be noted that this assessment of the relative 
importance of the various mechanisms is limited in application 
to the space and time scales associated with this model. 
However Ikeda et al . (1988) have presented similar conclusions 
based on a model of much smaller scale in the Labrador Sea. 
Additionally, Walsh et al. (1985) and Hibler and Bryan (1987) 
also present conclusions consistent with those here regarding 
the dominance of thermodynamic processes in the determination 
of ice cover extent. 

Correlations between simulated ice area and observed ice 
area were calculated for each major change in the model 
configuration. Two statistical approaches were used to 
determine if the difference between correlations obtained from 
each different model configuration were large enough to be 
statistically significant. Limited success was achieved in 
that the difference between the two model runs from the 
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interactive ocean experiment, in region 3, were 95% 
significant for both statistical approaches. However, 
statistical significance could not be determined at the 95% 
confidence level for the majority of regions because the 
variability of the correlations was too high. The large 
correlation variability indicated that the model was doing 
very well in some time periods but poorly in others. Two 
explanations were proposed. First, that some process 
controlling the ice area, for example the mixed layer, could 
be inadequately represented by the model. In those conditions 
where the process was relatively unimportant, the model would 
do fine. However, if conditions changed wherein the process 
became important, the model simulations would be degraded. 
A second possible explanation is that the atmospheric forcing 
is not being handled correctly. In particular, the monthly 
averaged wind stress is probably not accounting for the 
effects of short term, intense weather features. The dynamic 
mechanism experiment pointed out the large impact that the 
wind stress has on the ice field; therefore, this 
underestimation of the wind stress could result in periods of 
poor model performance. 

This work has shown that the relative importance of the 
various processes controlling the ice cover in the Arctic is 
regionally dependent. For example, in the western Arctic the 
dynamic effects assist the ocean heat flux in clearing the 
western boundary seas of ice. In contrast, the dynamics in 
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the eastern Arctic tend to oppose the clearing of ice from 
that region. That is, the dynamic forcing pushes the ice edge 
southward while the ocean heat flux tries to melt the MIZ 
northward. Various other combinations of effects are also 
evident in smaller sub-regions within the Arctic Basin. In 
order to allow for these variations, an ice model applied to 
the Arctic must be able to represent each mechanism reasonably 
well and be sensitive to the regional differences that change 
the importance of those mechanisms. Inclusion of a fully 
interactive prognostic ocean model in a linked ice-ocean model 
is necessary to do this. 

B. MODEL LIMITATIONS 

The discussion above already pointed out a potential 
problem if the monthly averaged wind fields did not provide 
sufficient temporal resolution of the wind to account for 
strong, episodic wind anomalies. The limitations of the 
forcing data must also be considered. This must include not 
only the atmospheric forcing, which is known to have some 
weaknesses, but also the prescribed inflow at the boundaries. 

A thermodynamic role for the ocean has been proposed 
wherein the ocean acts to precondition the waters below the 
surface layer in certain regions with heat. However, despite 
the recent reassessment of the Intermediate Water circulation 
by Aagaard (1988) which would appear to support the large 
scale circulation simulated by the model used here, 
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considerable ambiguity on the smaller scale remains. 
Furthermore, although the simulations indicate that 
considerable interannual variability exists for the 
geostrophic surface currents, the accuracy of this variability 
has not been determined. 

The representation of the mixed layer, which is the 
critical link between the ice and ocean portions of the model, 
was very simplified. It was noted above that the processes 
which occur above and below the mixed layer generally act to 
position the heat under areas of reduced ice thickness. 
However it is the action of the mixed layer, largely driven 
from the surface, which actually connects the ocean with the 
ice. The simulation of any such connection is currently 
limited by the vertical resolution of the upper ocean layers 
and the simple representation of the mixing which occurs 
there . 

The drag coefficients used here are the same for air/ice 
and air/water and constant for the ice/water. Recent studies 
in the Arctic have shown that these drag coefficients vary 
dramatically with the surface type, ice concentration, ridge 
concentration and several other factors (e.g., Guest and 
Davidson, 1987) . The drag coefficients determine the stress 
and are therefore important to the evolution of the ice field. 
Their simplification reduces the dynamic complexity of the 
forcing. 
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C. PROPOSED FUTURE WORK 

The surface mixed layer is the link between the ice and 
ocean sub-models, and consequently has an effect on the 
majority of the ice forcing parameters. The relative 
importance of many of these parameters has been examined in 
this work. The most active times of large-scale ice edge 
advance and retreat occur at the start of the freezing season, 
and the middle of the melt season, when the pack begins to 
break apart. It is also the same times when the mixed layer 
can undergo significant changes due to convective and dynamic 
mixing. Several authors have noted the importance of 
incorporating a realistic mixed layer model in model 
simulations of sea surface temperature (SST) . Fleming (1987) 
and Garcia (1988) have shown the strong correlation between 
SST and ice concentration in the Arctic. It would therefore 
be worthwhile to examine the importance of including a mixed 
layer model into the linked ice-ocean model used here. 

The role of atmospheric forcing in the interannual 
variability of Arctic sea ice was examined in Hibler and Walsh 
(1982). That paper concluded that simulations of interannual 
fluctuations in the ice cover were improved with the inclusion 
of interannually-varying, atmospherically-forced dynamics. 
A similar conclusion was drawn by Walsh et al., (1985) using 
a more elaborate model; however neither of these efforts 
included interannual variations of the ocean. A re-analysis 
of the role of interannually variable atmospheric forcing 
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would be valuable using this model, as it includes an 
interannually varying ocean. Further, mean monthly vice daily 
atmospheric forcing could be examined. 

The drag coefficients currently used are quite simple, as 
noted above. A sensitivity study examining the importance of 
the drag coefficients to simulations of the ice cover would 
be valuable and, if warranted, a more elaborate representation 
which accounts for the coefficient dependencies should be 
developed and incorporated into the model. It is anticipated 
that this would generate more leads in the ice, encouraging 
increased ice growth and a thicker ice cover. Greater 
compression might also occur in the MIZ and rough surface 
compression regions where the drag coefficients are larger. 

The work proposed above would progress our knowledge of 
ice cover simulation and Arctic Ocean circulation on the 
seasonal time scale. The next logical step would be 
incorporation of this model into an atmospheric GCM. This 
would be feasible using the current code because it has a 
relatively small computation requirement. The computation 
requirements of such a coupled GCM are expected to be well 
within the capabilities of the next generation of 
supercomputers . 

Finally, in order to improve the simulation and prediction 
of Arctic ice on a daily basis (a capability desired by most 
of the world's Navies), the effects of mesoscale ocean 
features and short duration atmospheric forcing must be 
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examined. To this end, the next major improvement of this 
model should include an increase in the horizontal resolution 
to at least 1/8 or 1/10 degree, bi-harmonic diffusion, an 
increase in the number of vertical levels to approximately 40 
and daily if not six hourly atmospheric forcing. 



D. FINAL GEM 

Walsh et al., (1985) stated: 

. . .model-derived trends (of Arctic ice cover) may be 
misleading in the absence of a realistic treatment of ice 
dynamics . 

To this I would add: 

and the direct and interactive effects of three- 
dimensional ocean circulation. 
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